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Preface 



The recent developments in dynamic cardiac imaging and modeling are the re- 
sult of an increasingly fruitful cooperation between theoreticians, engineers, and 
practitioners. The standard clinical parameters are progressively benefitting from 
more accurate quantitative measurements thanks to a simultaneous use of so- 
phisticated data acquisition techniques such as multichannel ECG and MEG, 
dynamic US, doppler, X-Ray, GT, helical GT, MRI, fMRI, SPEGT, PET... data 
and advanced physical and mathematical tools. It is becoming possible to effi- 
ciently combine prior anatomical and functional knowledge with advanced 3D 
spatio-temporal digital data. The challenge now is clearly to correlate the micro 
structure and function of the organ from the cellular level with its macroscopic 
functional behavior. 

The EIMH 2001 international workshop aimed to promote collaboration bet- 
ween scientists in signal and image processing, applied mathematics and physics, 
biomedical engineering and computer science, and experts in cardiology, radio- 
logy, biology, and physiology. The EIMH 2001 workshop, in its first year, focused 
on complex heart models involving anatomical and functional information. The 
goal is to gradually move toward hybrid 3D biomechanical and electrophysiolo- 
gical models able to simulate the dynamic behavior of the heart. 

We feel that EIMH 2001 presented an opportunity to create a cardio- vascular 
research network at a European level to federate the main pluridisciplinary 
groups involved in advanced research in cardiovascular imaging and modeling 
in Europe. The idea is to create what could be called the “European Beating 
Heart Project” . The purpose of such a project is to join research efforts to im- 
prove diagnosis and therapy of pathological dysfunctions of the heart, for the 
benefit of the patient, through fundamental and clinical research. 
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Abstract. Numerieal models of various aspeets of eardiae funetion are useful in 
simulating pathologies and therapies. Sueh models eannot be easily eombined 
beeause of differenees in input and output strueture. To harbor the different 
models we have developed a framework, eonsisting of a graph network with 
nodes (eompliant eavities) and edges (flow duets). Thus using a toolkit of 
modules faeilitates tailored eardiae modeling. 



1 Introduction 



Numerical models of cardiac mechanics are expected to be valuable tools in 
simulating cardiac pathologies and evaluating different strategies for therapy. In 
designing a cardiac simulation, one has to decide, which physiological mechanisms 
have to be included in a simulation. Important aspects may be modification of the 
depolarization wave by the degree of stretch of the myocardial tissue, coronary flow 
dynamics, the metabolic supply-demand balance of the myocardial tissue, the internal 
fiber structure within the heart, and adaptation of the structure to changes in load. 

In modeling cardiac function in more detail, interaction with the rest of the 
circulation is crucial. Focussing on the heart, a cardiac beat is initiated by a 
depolarization wave originating in sinus node, spreading out subsequently over the 
atria, the atrioventricular node, and both ventricles. As a consequence the heart 
contracts. Pressure rises in the ventricles. After exceeding pressure in the large 
arteries the outflow valves open. While the myocardium relaxes, flow diminishes and 
the outflow valves close. Then, the pressure head from the atria to the ventricular 
cavities causes the atrioventricular valve to open, resulting in filling. 

If regional myocardial function is a focus, the circulation is an important boundary 
condition to the mechanics of the wall. Venous pressures determine diastolic filling. 
Arterial pressures determine outflow phenomena. The balance of inflow and outflow 
determines the volumes in the cardiac cavities. When knowing the latter volumes, 
numerical models of cardiac mechanics may be used to determine global and local 
mechanical load of the myocardium as a function of time. 

Currently, most models of cardiac mechanics are focussed on one cavity, generally 
the left ventricle. Boundary conditions, set by the circulation to cavity mechanics, are 
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often simplified to a preload and afterload pressure level. One may think of 
simulating the mechanics of the various cavities with their mutual interactions during 
for instance pacing. These interactions are very sensitive to relative timing of events 
in the various cavities, thus needing to include the circulation in the model. 

Because most models are made with a specific goal, little attention has been paid 
on combining the strength of the various models. They are generally not suited to be 
mutually combined. On may solve this problem by simulating all cardiac and 
circulatory mechanisms in great detail and as accurate as possible according to current 
physiological knowledge. Despite optimism about increasing possibilities in 
computing, such approach is not likely to work. According to standard rules in 
modeling strategy, one should first determine the goal of the simulation. In reducing 
modeling efforts to a minimum a focus is chosen near which accuracy should be as 
high as required. Further from the focus accuracy may be less, thus enabling 
important savings of effort. A disadvantage of this approach is that each problem will 
result in a different model, generally not compatible with other models. 

The question arises how to combine the advantages of one comprehensive and 
accurate model with the tailored and much faster models. We propose to think in 
modules. The framework of the circulation model should be a working, simple and 
well-structured model of the circulation with the heart. When choosing a focus of 
simulation, more accurate ones should replace modules near the focus. When 
replacing a module, the structure of communication with the rest of the model should 
be left unaffected. Ideally, a toolkit of modules with standard interfacings should be 
available. For a specific problem, the appropriate tools should be selected and 
composed to a tailored, complete model system. An important advantage of this 
approach is that developed specific modules can be used in other applications and by 
other groups so that not much modeling work has to be repeated. 

In the present study we introduce a working model of the circulation, having a 
modular structure, and containing a very limited number of module types. The 
structure is that of a mathematical graph, containing nodes interconnected with edges. 
The nodes are compliant cavities, and the edges are flow ducts with inertia and 
resistance. A cardiac chamber module is a non-linear time dependent compliant 
cavity. A valve is a non-linear, flow-direction sensitive resistance with inertia. For 
each module the type is defined, the list of parameters required by the type function, 
and the interconnection with the other modules. Thus a lumped model of circulatory 
dynamics is obtained, leading to a set of differential equations. Cavities and ducts 
have as state variable volume and flow, respectively. The derivatives of the related 
state variables follow from the module functions. 

A few examples will be shown of focussing by replacing modules by more 
accurate ones, leaving the basic structure of the model intact. The examples concern 
the introduction of 1) sarcomere fiber dynamics, 2) transmural gradients in myofiber 
mechanics, and 3) finite element analysis of electrophysiological depolarization wave 
and ventricular mechanics. 



2 Setup of the Model 



The model of the circulation contains cavities as nodes, interconnected with flow 
ducts as edges, thus forming a graph network. The definition of the model is 
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transformed to a system of differential equations in a set of state variables. State 
variables are linked to cavities and ducts. Knowing the values of all state variables, 
the state of the system is fully determined. Solution of the system of equations 
requires that 1) the derivatives of the state variables are expressed as a function of the 
state variables themselves, and 2) the values of all state variables are known at the 
starting moment of the simulation. 



2.1 The Module of Type Cavity 

A node of the model represents a cavity V. Pressure p in this cavity is a function of its 
volume: 



p = function (v) ( 1 ) 

Most collagen based passive elastic biological cavities may be modeled by an 
exponential pressure volume relationship (7). In the circulatory system many cavities 
contain muscular tissues in the wall, among of which the cardiac chambers are the 
most obvious ones. The pressure p in such a chamber depends on cavity volume V 
and time t. The simplest model is a conventional passive elastic chamber, completed 
with a variable elastance model (8, 9): 



Parameters p , . and V ^ are linked to the passive element. Parameter 

A pas-^ pass,o pass a aci,u 

represents the small residual volume, at which no active pressure can be developed. 
The function E^^^(t) contains the time function of activation. The moment t=0 defines 
initiation of the contraction cycle. 

For the solution of the differential equation the time derivative of the cavity 
volume is expressed as the sum of flows q. in the ducts directed to the cavity: 

dV ^ (3) 



2.2 The Module of Type Duct 

An edge of the model represents a duct having inertia and resistance. Flow is a state 
variable. In modeling the duct, the time derivative of flow must be calculated as a 
function of flow and pressures in the proximal (Pp,^,,) and distal node (PdiJ. It holds: 

dq _ (4) 

function(ppj.Q^, Pdist’9) 

The pressures are a function of the volumes in the nodes by the pressure volume 
relation according to Eq. (1-2). Thus, for a duct the important condition is satisfied 
that the derivative of the state variable q is a function of state variables. Consider a 
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duct as an inertia in series with a resistive module. The pressure drop p^ over the non- 
linear resistive component is a function of flow q: 

Pj. = function (q) (5) 

If this function is linear, the resistive part of the flow duct is a linear resistance. The 
remaining pressure over the inertia is the pressure drop Pp^^^ - p^^^^ over the duct minus 
the resistive pressure drop p^. The inertia may be non-linear. For the flow derivative it 
holds: 



^ = function(pp,ox - Pdist “ Pr . q) 

Note that a bloodvessel segment may be modeled by a cascade of compliant 
cavities alternated with inertive ducts. Such model simulates realistic conduction 
properties of the arterial pulse wave. 

We have modeled a cardiac valve as a duct with inertia in series with a non-linear 
resistance. The physical background of the resistive part has been modeled as flow 
through an orifice with area A, which depends on current flow and on pressure drop. 
Furthermore, there may be a resistance R in series with the valve. The valve is 
considered open if either flow or pressure drop or both are positive. Thus Eq. (5) is 
specified to: 

_ P^'l^l I ^ A = ^ (Pprox ^ Pdist ) • ^open 

2A^ 1 else :Aig^k 

Note that Aj^^^ may be negligibly small, but not zero. The time derivative of flow is 
determined by the remaining pressure loss over the inertia, modeled by a flow channel 
with length 1 and area A. Thus Eq. (6) is specified to: 

dq _ A(pprox- Pdist -Pr) (8) 

dt pi 

Parameters p, R, A^p^^, and 1 are parameters to be provided. 

If a duct is a resistance without inertia, flow itself is not a state variable in the set 
of differential equations, but it is a straightforward function of the state variables of 
the proximal and distal cavities: 

^ “ (pprox “ Pdist )/^duct 



2.3 Implementation 

The structure of cavities interconnected with ducts is convenient as long as cavity 
volumes are state variables. It is inconvenient to have nodes not having volume as a 
state variable because then duct flows become mutually inconveniently coupled. To 
circumvent such problems it may useful to attribute some compliance to the cavity. 
The related time constants in the differential equation should be kept sufficiently 
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small to keep phase errors negligible and sufficiently large to maintain numerical 
stability. 

When setting up the simulation, it is important to decide what minimum rise time 
will be allowed for the dynamical behavior. For numerical stability all modules 
should have time constants not shorter than this rise time. 



2.4 Design of the Circulation Model 

The dynamic model of the normal circulation contains a systemic and a pulmonary 
tract (fig 1). Ducts interconnect cavities. This model will be used as the basic 
structure of cardiovascular modeling. Specific features can be focussed on by 
replacing modules by more detailed ones, while maintaining the original links. New 
ducts can be added to simulate many types of shunts, such as fistulas and cardiac 
shunts in newborns. 



Left heart 



Systemic Circulation 



sa sRpvvw sv svRWfr- 




Right heart 



Pulmonary Circulation 



Fig. 1. The circulation has been modeled by a left and a right circuit of compliant cavities, 
interconnected by ducts, such as valves and resistances. The dashed arrow indicates a shortcut, 
used when focussing on simulating left ventricular hemodynamics. Then the LV module was 
replaced by a more accurate one. 3 Specific Simulations 



3 Specific Simulations 

The model setup has been used to simulate the boundary conditions set to models, 
focussing on cardiac mechanics. As examples we have chosen 1) replacing the 
variable elastance model by a first order stress-strain model of the sarcomere, and 
coupling of sarcomere dynamics to left ventricular cavity mechanics, 2) relating 
cardiac deformation as measured with MRI with regional myofiber mechanics, 3) the 
effects of changing the pacing site on the ventricle for the sequence of local myofiber 
contraction. 



3.1 Linking Ventricular Dynamics to Sarcomere Dynamics 

In the model of Fig. 1, pressure of the left ventricle is calculated from its volume by 
an equation of type Eq. (1). In linking ventricular dynamics to sarcomere dynamics, 
local differences are considered not to exist. So we used a simple 1 -fiber model as 
presented earlier (1). The strength of this model is its accurate description of the non- 
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linear relationships between cavity volume and pressure p^^^ on the one hand, and 
sarcomere length 1^ and myofiber stress s^ on the other. It holds by good 
approximation: 



ls=l 



sO 



1 + 3- 



V., 



X 1/ 

^/3 



V. 



wall 



and Pcav = 



1 + 3- 



V, 



V. 



wall 



Sf(ls.t) 



( 10 ) 



Parameter represents wall volume and 1^^ is a scaling length for sarcomere length 
in the wall. In maintaining the links with the circulation model, Eq. (10) is used to 
express pressure p^^^ as a function of cavity volume given that the expression 
s/l^,t) for myofiber stress is known. For the latter expression we have used a first 
order Maxwell model, consisting of a passive elastic element in parallel with a series 
combination of a contractile element with length l^j, and a series elastic element l^j,. 
Velocity of contraction dl^.j,/dt for the contractile element depends on external load 
and on the degree of activation, incorporated in the function Sj^^(t). We have used: 

^pass “ ^passO ij — Sjg^ (t)(lg — l(^g )/lggQ ^ ^ 

~ ^pass ^act ~ ~ (^s ~ ^CE )/^SE0 ) ^max 

Parameters k, 1^^, and function s.^^(t) were derived from physiological 

experiments on isolated cardiac muscle. From Eq. (11) the length of the contractile 
element appears an additional state variable, which is solved in the framework of the 
whole circulation. In Fig. 2 the results of a simulation are shown for occurring 
pressures and flows in the heart and circulation. Since the hemodynamics of the left 
ventricle was in focus, a shortcut was made for blood flow from the systemic veins to 
the mitral valve, thus excluding the pulmonary circulation (Fig. 1). 




Fig. 2. Simulation of normal hemodynamic variables, using a 1 -fiber model of left ventricular 
mechanics, applying Eq. (10-11). 



3.2 Analysis of Cardiac Deformation as Measured with MRI 

A method was developed to assess transmural differences in cardiac function. With 
MRI-tagging motion of the heart can be measured in a number of slices. When 
imaging two parallel short axis cross-sections near the equator, the most obvious 
modes of deformation are contraction and torsion. Torsion is defined as a base to apex 
gradient in rotation angle of the cross-section of the left ventricle around the base to 
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apex axis. Because of the helical structure of the myofibers in the cardiac wall, this 
torsion contains information about the transmural course of myofiber shortening. The 
simplest way to model these effects is considering part of the left ventricular wall to 
be a segment of a cylinder (2, 4). In the model of the whole circulation, the ventricular 
model can be implemented as follows (Fig. 3). The deformation of the cylinder is a 
function of left ventricular volume Vj^, torsion T, and base to apex length L. Given the 
length of the contractile elements in the sarcomeres within the wall, stress was 
calculated (see also Eq. 11). Equilibria of forces are used to solve T and L. Thus, 
geometry and sarcomere length inside the wall are known. Left ventricular pressure is 
calculated from the related stresses. Besides, these stresses also determine the time 
derivatives of contractile element lengths in the wall, which are the state variables 
associated with the myofibers in the wall. Thus, pressure is a function of volume, and 
the derivatives of the introduced state variables are known. 




Fig. 3. Pressure as a funetion of eavity volume. Length and torsion determine sareomere length 
L From the resulting stress follows a foree balanee, to be satisfied by proper ehoiee of torsion 
and length. Cavity pressure results. Contraetile element length l^^is a loeal state variable, 
needing its loeal derivative. 



3.3 Sequence of Mechanical Activation During Pacing 

Normally, a heartbeat initiated by depolarization of the sinoatrial node (SA). 
Subsequently, an electrical depolarization wave crosses the atria. With some delay, 
the wave travels through a narrow duct, the atrioventricular node, thus reaching the 
ventricular conduction system. Via the bundles of His and the Purkinje conduction 
system, the ventricles are depolarized. Depolarization of the cardiac cell initiates 
contraction of the myofibers. So, the sequence of electrical depolarization is reflected 
in the sequence of mechanical contraction. One of the strategies to obtain non- 
invasively the sequence of cardiac activation is the measurement of cardiac 
deformation, followed by a reconstruction. For this reconstruction, a model of cardiac 
mechanics is used. 

To simulate the sequence of electrical (6) and mechanical activation, the heart 
should be modeled in 3 dimensions. Conventionally the deformation of an object is 
solved numerically by the finite element method (FEM) (3, 5). With this technique, 
forces are exerted to the nodes of the mesh. Deformation of the heart is adapted 
numerically by iteration until the equilibrium of forces is satisfied. So, for a known 
pressure load, geometry and cavity volume are calculated. Unfortunately, this 
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condition is not compatible with our model setup for the circulation, because pressure 
should be calculated as a function of volume. A regular solution to this problem is 
iteration of pressure until FEM-volume fits the volume in the circulation 
model (Fig. 4). Such an approach is however very time consuming, because the FEM 
is incorporated in the iteration. Alternatively, one may use the 1 -fiber model (Eq. 10- 
1 1) to accurately predict left ventricular pressure. The difference between and 
can be used to adjust the parameters of the 1 -fiber model to get still better estimates 
for left ventricular pressure. Correction parameters may be implemented as state 
variables to be integrated in time. 



Fig. 4. Implementation of a FEM model of eardiae meehanies in the eireulation model, eavity 
pressure should be deseribed as a funetion of eavity volume Vj^. The FEM model generates 
volume as a funetion of pressure. The problem ean be solved, either by iteration (left) or using a 
1 -fiber model, replaeing the FEM, while adjusting the 1 fiber model by feedbaek. 
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4 Conclusions 

A model of the circulation can be used conveniently to harbor modules describing 
heart mechanics. The structure of the circulation model is a graph network with 
cavities as nodes and flow ducts as edges. The cavities contain volume as a state 
variable for the governing set of differential equations. In focussing on various 
aspects of heart or circulatory dynamics, tailored modules can be inserted, while 
maintaining overall structure of the model. A toolkit of modules may be extended 
flexibly. 
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Abstract. This paper introduees a global geometrieal modeling approaeh to rep- 
resent the heart and its main vessels. It takes into aeeount that this model will be 
used for heart motion analysis: we have emphasized the role of organ struetures 
beeause they should be a support for motion representation. But heart and ves- 
sels refer to different types of struetures: thus, we have defined two different 
geometrieal modeling approaehes to handle these two anatomieal elements. 

Keywords. 3D medieal imaging, heart motion modeling, reeonstruetion 



1 Introduction 

Geometrical modeling of the heart and its main vessels is required as a support for 
analyzing their motion. In this paper, we suppose that a segmentation process has been 
performed and has produced a binary volume, our goal being to associate a geometri- 
cal model to it. 

An organ can be fully represented by the closed surface that bounds it. This repre- 
sentation is more synthetic than a volumetric description, but we may want to know 
easily if a given point is located inside or outside the organ: the geometrical model 
must facilitate the answer to this question. 

The geometrical model must also give indications on the organ’s topology and 
morphology, and especially on its structure, such as on its local variation properties. 

The modeling approaches we propose to reconstruct the heart and its main vessels 
satisfy these constraints. 



2 Heart Modeling 

2.1 Our First Approach 

Our first goal was to characterize the global shape of the heart through its structure. 
For this reason, we have decided to orient our research work towards models based on 
implicit surfaces, their skeleton giving interesting information on their morphology. 
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The reconstruction process is based on the use of the medial axis [1] as a potential 
source. A specific potential is attached to each point of the medial axis so that the sum 
of these potentials provides equipotential surfaces [2,3], and especially one surface 
which has to represent the binary volume’s crust as well as possible [4]. 

The potential function [5] depends on three parameters r, R and k that completely 
characterize it. Its description is given on Fig. 1. 
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Fig. 1. The potential function and its formal definition 
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Fig. 2. The digital volume of a fetus’ heart, the medial axis, its reduction and the reconstruction 
by implicit surfaces 

The reconstruction process consists in determining the set of parameters that mini- 
mizes a given energetic criterion. This criterion characterizes the similarity between 
the equipotential surface and the crust of the binary volume. The entire process is 
described in [6]. Fig. 2 shows an example of a heart reconstruction based on this ap- 
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proach. The delay for the reconstruction process is about 10 minutes on a 450MHz 
Pentium III processor (512 Mo RAM). It seems to be satisfactory within a medical 
examination. 

Such a geometrical model provides a surface and volume description (the volume is 
obtained by the value of the global potential). It also provides an interesting data com- 
pression, but two problems rise from this approach: 

• the number of potential sources is too high, even after medial axis reduction 

• there is no structure underlying the set of potential sources 

We found a solution to these two problems through a multilevel approach which is 
based on the same principles but takes into account strong structural elements. 



2.2 A Multilevel Approach 



In order to find the binary volume’s structure, we have analyzed this data set through 
different levels of resolution. The first level, which is very rough, provides the global 
structure of the shape. In fact, the result is made of a few points and it is easy to de- 
termine the structure of such a data set (although it is not explicitly used in the recon- 
struction process). This set of "inner" potential sources gives a basic volume that will 
be covered and enriched by refining the resolution and completing the external layer. 
Fig. 3 illustrates this process on a simple example. Each case refers to a sub-resolution 
of the starting digital volume, strictly included in the latter. The object is seen as a set 
of layers, in order of structural importance. 






Fig. 3. The layers of an objeef s model 




This multilevel approach is very powerful for many reasons: 

• it provides an efficient data compression; 

• it gives a representation of both surface and volume; 

• it enables the understanding of the heart structure through different layers and 
skeletons (each one of them is related to a specific layer); 
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• it provides an interesting support for motion analysis because motion can be split 
up into complementary components (global motion associated with the inner 
structure, etc.). 



3 Vascular Structure Modeling 

The approach described in the previous section cannot be efficiently used for vascular 
structure modeling because it does not reflect an over structure related to the global 
orientation of vessels (tree-like structure). Thus, we have developed a complementary 
model to represent such specific anatomical structures. 



3.1 Characterization of a Tree-Like Free-Form Surface 

Most of the research work done in this field describes the global surface using a set of 
patches with conditions on their boundary so as to obtain a given level of continuity 
globally [7,8,9,10,11]. 

We have chosen a different approach that consists in using a single parametric 
patch on which we set topological and morphological conditions. 



3.2 Use of a Parametric Model 

The shape of a surface is characterized through its topological, morphological and 
differential properties. Using a classical parametric surface as a basis for the modeling 
approach enables the control of its differential properties (e.g. using a cubic B-Spline 
surface produces a model). 

Topological properties are controlled through the definition of invalidation areas in 
the parametric domain: for example, a bounded parametric domain with invalidation 
areas made of two connected components (i.e. two "holes") is homeomorph to a 
sphere with three holes. Thus, a tree-like surface with N ending parts is topologically 
equivalent to a single patch with (N-l) holes, or to a limited parametric domain with 
(N-l) invalidation areas. 

The main difficulty in using such an approach is to take into account the morpho- 
logical properties of the surface. And this is the object of our developments. As an 
example, we show the drawing of an isoparametric line on a surface which is a single 
patch without hole but with an underlying structure and we note that isoparametric 
lines do not reflect this structure (Fig. 4a). Even on such a very simple example, we 
would like these isoparametric lines to be associated with sections (Fig. 4b). 

We can provide such a morphological control on the surface through an intermedi- 
ate interface domain level that reflects the tree-like structure through its parameteriza- 
tion. 
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Fig. 4. a. Image on the surfaee of an isoparametrie line of the definition domain 
b. A eoherent parameterization of the domain aeeording to the surfaee 



3.3 Need for an Intermediate Interface Level 

The final surface can come from any model which uses control points and a classical 
parameterization. But we do not wish to access directly to these control points to ma- 
nipulate the surface. Instead, we define an intermediate level of domain with another 
parameterization. Then, the surface defined by this intermediate level is used to sam- 
ple the final surface’s control points. 

Let us consider a junction (the simplest tree-like structure) as an illustration of our 
approach. The parametrical domain has two invalidation areas that are obtained 
through the definition of potential sources (we invalidate the areas where the global 
potential is over a given value, or under another given value). Equipotential lines can 
then be associated with sections of the junction and provide a new parametric ap- 
proach (one parameter is defined along these lines and the other in their orthogonal 
direction), as illustrated on Fig. 5. 



a 




Fig. 5. a. The parametric domain of a branching shape 

b. Isoparametric grid of a branching shape domain 



New control points are then obtained through the new parameterization of the in- 
termediate domain. Fig. 6a shows the correspondence between these two levels of 
control points. We have illustrated our approach through a junction but it is easily 
generalized to any tree-like structure by associating the basic structure to imbedded 
potential sources (Fig. 6b shows the structure and Fig. 6c shows the potential value 
that is used to produce the new parametric domain). Fig. 6d illustrates the use of the 
model for a vascular structure. 
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Fig. 6. a. Computation of the control points from the meta control points grid 

b. Tree-like free-form surface of a branching shape 

c. Parametric domain related to the structure in b 

d. Example of a vascular structure 



4 Results and Future Work 

These two approaches are complementary and well adapted to each case (heart and 
vessels). Because vessels are not surfaces, we need to define their two walls as sur- 
faces: the geometrical modeling approach described earlier enables such a representa- 
tion by using two of these surfaces and a correspondence between them (the structure 
is the same, as with the parametric field and control point definition, the only differ- 
ence is the 3D control point location), this correspondence facilitating the description 
of the vessel tissue’s thickness. 




Dtghal volume Implicit surface (Heart] 



Fig. 7. Heart’s data, its reconstruction with implicit surfaces and aorta model 
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These two geometrical modeling approaches are also interesting for motion de- 
scription because we can try to analyze and to encode this information at different 
levels (structural level, surface level). This could be done by making a link with a 
physical animated model of the heart. Thus, one aim would be to compute functional 
parameters such as the myocardial mass or the volume variation through time. 

One of the remaining problems - and it is a main research axis for us - is the study 
of the global model’s coherence, especially at their junction (aorta, for example) or 
when they are in contact (coronaries, for example). 
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Abstract. Individual 3D boundary element models ean be used in solving 
inverse problems in eleetro- and magnetoeardiographie measurements. In some 
eases 3D data, sueh as Magnetie Resonanee (MR) or Computed Tomography 
(CT) images, are not available. Therefore, it would be useful to be able to use 
2D images sueh as X-ray projeetions for ereating 3D models. The aim of this 
work was to develop a software paekage for ereating a 3D boundary element 
heart model from two orthogonal X-ray projeetions. The biplane fluoroseopy 
images from a patient are digitized and the images are enhaneed with different 
image proeessing teehniques. The patient heart outline is segmented from the 
X-ray projeetions. The outline is eompared with virtual X-ray projeetions 
ereated from a prior 3D model segmented from MR images. The differenee 
between the outlines is used to deform the prior model. The quality of the 
digitized X-ray projeetions was notieeably improved and thus the heart outline 
segmentation was faeilitated. The deformation method implemented is robust 
and provides good results even when the souree parameters eontain errors. 



1 Introduction 

Individual 3D boundary element models can be used in solving inverse problems in 
electro- and magnetoeardiographie measurements. Many other applications exist for 
individual 3D models in the medical field. In some cases, Magnetic Resonance (MR) 
or Computed Tomography (CT) images are not available or the cost of obtaining them 
is too high. Modalities producing 2D data, such as X-ray projections and ultrasound 
images are more readily available. Therefore, it would be useful to be able to use 
these images for creating 3D models. 

One or two X-ray projections of an object are usually not enough to reconstruct a 
3D object. Several X-ray projections of the same object from different angles are 
needed to create a good 3D approximation of the object. In CT, hundreds of 
projections are used. If only a few X-ray projections are used, some a priori 
knowledge can be incorporated to improve the reconstruction accuracy. For human 
anatomy, the prior knowledge, i.e. a representation of typical human anatomy, can be 
extracted, for example, from MR volumes. This work is based on a method where a 
3D model was reconstructed from two 2D silhouettes using a prior 3D model [1]. 
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The aim of this work was to develop an easy to use software package for creating a 
3D heart model from two orthogonal X-ray projections. In addition, 3D model 
coordinates can be displayed on the original X-ray projections in 2D for persons who 
are used to locating objects from X-ray projections. Also, 2D coordinates, such as 
catheter locations, can be visualized with the 3D model. 

The developed application is not limited to creating only heart models but it can be 
applied more generally. Once a prior model of an object and two images of a new 
slightly different instance of the object are available, this application can be used to 
reconstruct the 3D geometry of the object. 



2 Methods 

The geometric prior model used in this work was built from MR images. The data 
volume used here originates from a cardiac cine MR sequence at diastolic phase. The 
geometric model is built by first segmenting [2] and then triangulating [3] the cine 
MR volume. The segmentation is based on free-form deformation of a geometric and 
topologic prior model. The model is matched to the edges present in the data volume. 
After segmentation, a set of points is selected from the segmented volume and the 
points are linked with each other in such a way that triangles are produced. The 
Voronoi-Delaunay duality is used in the triangulation. The epicardium and both 
ventricles are included in the model. 

The source images are sequences of X-ray images of the heart of a patient, taken at 
two different orthogonal angles. The images used in this work are taken during 
catheter operations and contain, in addition to the anatomical features, the electrodes 
and catheters used during the operation. The images are converted from analog video 
data into digital format with an image grabbing system. 

In order to improve the low quality of the source images, a number of image- 
processing algorithms were implemented. The source image frames are combined 
with intelligent averaging to reduce the image noise and interleaving to enhance the 
resolution. A specifically modified local contrast enhancement algorithm is applied, 
where the nature (for example, the constant black border) of the source images is 
taken into account. Local image neighborhood statistics are compared with global 
statistics and pixel intensities are modified to increase local contrast. The 
improvement in image quality helps the segmentation and visual inspection of the 
image. An example of the image enhancement is displayed in (Fig. 1). 

Several methods have been proposed for automatic segmentation of X-ray 
projections, such as active contour models [4] and united snakes [5]. The quality and 
type of the source images, used in this work, makes them very difficult for automatic 
segmentation. The problem is further complicated by the fact that the object is often 
seen only partially in the images. A method, introduced in [1], tries to overcome 
these problems using prior model information and elastic matching. However, this 
method has not yet been tried for the type of source images used in this work. The 
final segmentation of the heart outline is done manually by looking at the pre- 
processed source image. 
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Fig. 1. Image enhancement: (a) Original digitized fluoroscopy image; (b) Image processed with 
frame interleaving, image averaging and local contrast enhancement. The anatomic features in 
the image and the catheters are now more clearly visible. 



The outline is compared with virtual X-ray projections created from a prior 3D 
model, carefully segmented from MR images of the heart. The 2D vectors generated 
by matching the actual outlines with the prior model outlines are back-projected onto 
the 3D model surface. Vectors are interpolated for each model surface point. The 
original model points are moved according to the displacement vectors. The 
ventricles of the heart, not visible in the X-ray projections used, are deformed using 
local scaled averages of the deformations calculated for the epicardium (Fig. 2). 

As a result of the previous stages, necessary information exists to map 3D coordinates 
back to the original 2D projection coordinates. This information can be used to 
display 3D points on the original 2D projections, for example, on the ventricle 
surfaces, such as Josephson points. The points can be color-coded to display 
additional information. Also, 2D coordinates, such as catheter tip locations, on the 
original projections can be visualized in 3D with the new model (Fig. 3). 



3 Results 

The quality of the digitized X-ray projections was noticeably improved and thus the 
heart outline segmentation was facilitated (Fig. 1). At this point the relevant features 
of the image are considerably more prominent than in the original X-ray image. 

The direct validation of the results with actual cases where both 2D silhouette 
images of the heart and otherwise generated 3D images of the heart could not be done 
since no such data was available. However, the results can be indirectly observed 
with just 3D data available. Silhouette images can be generated from the 3D data in 
various ways, for example, with virtual X-ray projections. These 2D images can then 
be used to generate 3D models with this application and the deformed models can 
then be compared with the original models. Such a comparison has been carried out 
for torso models in [1], where the mean error for the torso models generated was 
found to be approximately 1.2 voxels for a 128 x 128 x 100 volume. However, the 
geometric accuracy of the model is not very critical in this work, since the deformed 
3D model is inspected only visually. 
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Fig. 2. The deforming proeess: (a) The original heart model displayed from two orthogonal 
angles; (b) The patient X-rays where the new heart outline is segmented from. The outlines are 
eompared with virtual X-rays of the original model; (e) The 2D veetors generated from elastie 
matehing between the aetual outlines and the model outlines are moved onto the 3D model 
surfaee; (d) Veetors are interpolated for eaeh model surfaee point; (e) The new deformed 
model. 
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Fig. 3. Two orthogonal views of catheter tip locations digitized from 2D fluoroscopy images 
visualized in 3D with a heart model including the heart ventricles. 

The direct validation of the results with actual cases where both 2D silhouette 
images of the heart and otherwise generated 3D images of the heart could not be done 
since no such data was available. However, the results can be indirectly observed 
with just 3D data available. Silhouette images can be generated from the 3D data in 
various ways, for example, with virtual X-ray projections. These 2D images can then 
be used to generate 3D models with this application and the deformed models can 
then be compared with the original models. Such a comparison has been carried out 
for torso models in [1], where the mean error for the torso models generated was 
found to be approximately 1.2 voxels for a 128 x 128 x 100 volume. However, the 
geometric accuracy of the model is not very critical in this work, since the deformed 
3D model is inspected only visually. 

The accuracy of the results could be defined using similar simulations to the paper 
[1], as mentioned above. Because no information is available about the ventricles in 
X-ray images used, the interpretation of the ventricle surfaces is, however, only 
qualitative. Therefore, the geometric accuracy of the ventricles is not emphasized and 
studied in this work. Nevertheless, it is worth noting that since the prior model is 
used, the geometric relations between the epicardium and ventricles remain 
approximately correct during the deformation (Fig. 3). If a contrast agent is applied 
during the catheter operation, so that the ventricle becomes visible in the X-ray 
projections, a more accurate model of the ventricle can be reconstructed with the 
software. 

The information on the source imaging angle is manually entered by the user. 
Changes in the resulting output when the source angles are erroneously entered are 
shown in Table 1. A model was created using accurate orthogonal angles and it was 
set as the base for comparisons. The source angles for one or both views were 
modified and the resulting model was compared to the base model. The difference in 
entered and real angles and the difference between the result and correct object 
volume is displayed in the table along with the overlap for the result and correct 
volumes. When the patient heart corresponds well to the used original heart model, 
the method is very robust to errors in determining the source image angles. Since the 
goal of the model creation is to reconstruct the correct model shape and not the size, 
the actual simulation parameters need not closely correspond to the actual X-ray 
imaging parameters used in the source X-ray projections. The generated model can 
be easily scaled using, for example, information from ultrasound images. 
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Table 1. Changes in produced 3D model volume with different source imaging angles 



Angle difference 


Volume difference 


Volume overlap 


U 


+0.043 % 


99.2 % 


2° 


+0,061 % 


99.2 % 


10° 


-0,034 % 


98.2 % 


o 

O 

(N 


+0,040 % 


97.8 % 



The deformation method implemented is robust and provides good results even 
when the source parameters contain errors. The application is easy to use and new 
models can be created on-line with a common desktop PC. Currently the software is 
used to create individual anatomic heart models. However, the software is not limited 
to creating just heart models. 



4 Discussion 

Presently, only one object outline can be displayed and edited at one time. If multiple 
object outlines could be edited at the same time, this would simplify and speed up 
handling of multiple objects. Also, the relations between the objects could be more 
easily viewed and edited. 

Manual segmentation works well in practice and it can be accomplished in a few 
minutes. Automatic segmentation such as outlined in [1] could be implemented and 
tested. However, the visibility of the edges of interest in our images is much poorer 
than in the thorax- images used in [1]. Moreover, strong edges produced by catheters 
may distort automatic segmentation. A form of snakes [5] was used to segment the 
images but this relays heavily on the initial guess for the heart outline to be quite 
close to the actual outline for acceptable results. Snakes can however be used to 
segment the heart outline at different phases semi-automatically once the initial 
segmentation from one frame is accomplished. 

Only X-ray images were used as source images but other imaging modalities such 
as ultra sound imaging where the object silhouette can be seen could be used. Any 
2D image where the object silhouette can be observed could be used as a source 
image. Possible changes in the imaging geometry should be taken into consideration. 

An option for increasing the number of source projections could be implemented. 
In cases, where more projections are available, a more accurate model could be 
generated using all available data. 
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Abstract. In echocardiography, the radio- frequency (RF) image is a 
rich source of information about the investigated tissues. Nevertheless, 
very few works are dedicated to boundary detection based on the RF 
image, as opposed to envelope image. In this paper, we investigate the 
feasibility and limitations of boundary detection in echocardiographic 
images based on the spectral contents of the RF signal. Using the sys- 
tem approach, we study on models and simulations how the spectral 
contents can be used for boundary detection. We then introduce an orig- 
inal method of spectral estimation for boundary detection, and several 
images are analyzed with its mean. It is shown that, under the condition 
of high acquisition frequency, it is possible to use the spectral contents 
for boundary detection, and that improvement can be expected with re- 
spect to traditional methods. The conclusions may enable development 
of a robust boundary detection method, based both on the envelope and 
the spectral contents of the RF signal. 



1 Introduction 

The interest of echography is well known. As compared to other medical imaging 
methods it is a relatively low-cost, fast and non-invasive method. Applied to 
the cardiac field, it enables real time imaging of the heart. In the perspective 
of accurate assessment of the cardiac function, automatic boundary detection in 
echographic images appears as a key issue. This problem has been widely treated 
in literature, however, due to speckle noise and image artifacts, it still remains 
a challenging and open problem. 

The unprocessed radio- frequency (RF) signal from the ultrasound scanner is 
potentially rich in information about the investigated organ and the boundaries 
between tissues. Recent works in echography successfully use spectrum analy- 
sis of RF signals to characterize the explored tissues [1] [2]. However, none of 
those works apply spectrum analysis to boundary detection, which continues to 
be performed almost exclusively on envelope (B-scan) images (see for example 
[3]). Such images are obtained by demodulating the RF signal, and its spectral 
contents is lost in this process. 

In the current paper, we address the feasibility of boundary detection based 
on the spectral contents of the RF signal. With a Discontinuity Adaptive method 
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of spectral estimation introduced by us in [ 6 ], we investigate the advantages 
and limitations of such an approach, as compared to traditional methods based 
exclusively on the envelope of the signal. This work is carried out in view of its 
application to the field of cardiac imaging. 

2 Boundary Detection Method 

Due to the stochastic nature of the RF signal, exploring its spectral contents 
requires a spectral estimation method, and consequently, choosing appropriate 
spectral parameters: discrepancy between such parameters issued from different 
tissues will enable correct boundary detection. It has been proposed in [4] to 
model portions of the RF signal as an autoregressive (AR) process of order p, 
characterized by the variance of the driving noise and the reflection coeffi- 
cients 7 /e. Such a description is equivalent to the Power Spectrum Density (PSD) 
[5]. We have introduced in [ 6 ] an original method, combining spectral AR (para- 
metric) estimation with Discontinuity Adaptive (DA) smoothing. Other works 
in this area [4] [12] do not take into account discontinuities while smoothing, and 
therefore are less suitable for the task of recovering boundaries. Our method can 
be directly applied to boundary detection, based both on the envelope and on 
the spectral contents of the echographic RF image. We briefly recall the method 
here, and apply it in section 5 to detect boundaries on several echographic im- 
ages. 

The RF signal along each acquisition line is partitioned into short windows. 
At each consecutive AR order k we adopt the following two-step algorithm: (I) 
the reflection coefficients of order k {'jk) are estimated for all the windows of the 
image using the Burg iterative scheme, (II) based on the field iA of coefficients 
7 /c, the smoothed field iA is obtained by means of a DA method. This second step 
is essential in the context of boundary detection. The solution iA is obtained by 
minimizing a cost function of the form: 

C{fk I A) = ^ a, II 7* - 7| ||2 +r] ^ ^(7! - 7* ) ( 1 ) 

Vs s,tev 

The first term stands for data fidelity, where we consider the field iA as the 
observation data. The second term takes into account the smoothness constraint, 
formalized through Markovian Random Fields applied on a neighborhood V. 
A non quadratic ^-function results in smoothing relatively even areas of the 
image, while preserving discontinuities [7] [ 8 ]. The parameter 7 is a weighting 
term which tunes the conditions of smoothing. We use an equivalent description 
of the ^-function, based on the so-called dual fields [ 8 ]. We thus obtain a dual 
field associated with the smoothed field of coefficients, and it is this very dual 
field that we regard as a map of discontinuities. 

3 Theoretical Modeling 

We model the RF signal formation process by the system approach, largely used 
in literature [10] [11] [9]. The PSD of the RF image |/(/,p)p at the observed 
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point p is obtained by the multiplication of two terms: the pulse response at the 
focal point and the tissue response |T(/,p)p: 

\I{f,p)\^ = \H{f)\^-\T{f,p)\^ (2) 

Backscattering of the tissue is proportional to the 2^^ derivative of its acoustic 
impedance fluctuations in the direction of the wave propagation. In the adopted 
inhomogeneous continuum approach, the tissue is modeled as an ensemble of 
gaussian scatterers, randomly distributed in the tissue. The pulse, H{x,y), is 
described as a cosine wave of central frequency /o modulated by a gaussian. 

For the purpose of our simulations we define three tissues, determined by the 
shapes of the scatterers: (A) symmetric gaussian shape, of equivalent diameter 
d = 150/rm, (B) gaussian shape oriented at an angle (j) = 10^ of equivalent length 
di = 100/rm and diameter dt = lOjum and (C) very small gaussian scatterers. 
The corresponding tissue responses are given below, where k is the wave number 
and the term stands for the second derivative 

| T ^(/ c )|2 = ( 3 ) 

\TB{k)'\^ = . g-2fe^(d? cos^ 0+d?sin^ 0)/4 

Using the same framework, 2D acquisitions can be simulated [11]: 2D tissues 
corresponding to ^4-Care depicted in Fig. 1. These simulations mimic 3 tissues 
encountered in cardiac imaging: the external tissue, the cardiac muscle with its 
oriented fibers and blood with small blood cells. The task of boundary detection 
in cardiac imaging requires discriminating between such tissues. We simulate two 
acquisition pulses whose central frequencies /o of 2 and 5 MHz approximately 
correspond to the lower and upper bounds in human echocardiography: 

|P(/)|2 (5) 

The bandwidth of the pulse, introduced through the parameter cTj, increases 
with frequency. 




Fig. 1. (a) - (c) modeled tissues, corresponding to A-C (see text), (d) a scheme speci- 
fying the shape of the simulated body and the positions of the tissues 
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4 Modeling Results 



Based on considerations from section 2 , theoretical expressions of PSD of the 
resulting RF signals are derived. Fig. 2(a), (b) present the PSD of the signals, 
each one issued from insonification of one tissue (tissues 4, R, (7) with the pulse 
of 2 MHz and of 5 MHz. The three spectra are almost identical at acquisition 
frequency 2 MHz, and much more dissimilar at 5 MHz. Hence, it seems that 
segmentation will be much easier in the latter case. 

We subsequently simulate echographic imaging of the synthetic body (Fig. 1 ) 
with the two aforementioned pulses, obtaining thus two RF images. In regions 
corresponding to each of the tissues we estimate one spectral parameter: the 
reflection coefficient of order 1 ( 71 ). It is a complex coefficient, whose phase is 
dependant on the central frequency of the PSD. The higher the central frequency 
of the PSD, the higher the phase of 71 . Let us also remind that the phase 
rises in the anti-clockwise direction in complex coordinates. Due to the non- 
stationarity of the signal, estimation is performed on short windows, yielding 
several estimates of 71 . The results are represented in Fig. 2(c), (d) in complex 
coordinates. The dispersion of the values inside one population results from the 
stochastic nature of the tissue simulation [ 11 ], as well as from the process of 
numerical estimation itself. 




rc 








K 



(c) 



(d) 



Fig. 2. (a), (b): theoretical spectra (in dB) of RF signals acquired at /o, respectively, 2 
and 5 MHz. (c), (d) 3 populations of coefficient 71 represented in complex coordinates, 
with complex circle of radius 1, acquired at /o 2 and 5 MHz 
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It is clear that the three populations of reflection coefficients are far better 
distinguishable with an acquisition at higher frequency. This confirms previous 
observations. Detection of boundaries based on the spectral parameter 71 ap- 
pears feasible for the higher acquisition frequency, while it seems an extremely 
difficult task for the 2-MHz acquisition. 

5 Boundary Detection Results 

In order to confirm the observations of section 4, we process echographic images 
acquired in different conditions: a simulation, a human and a canine heart. For 
each case, we present several figures: the envelope image; boundaries detected 
with our DA method based on the envelope image, which corresponds to pro- 
cessing a conventional B-scan image; boundaries detected with our DA method 
based on the spectral parameter 71 and based on the spectral parameter 72 . 
Running our DA-method requires adjusting two hyperparameters: the weighting 
parameter 7 (see eq. 1) and a threshold parameter 6 for the function For each 
image separately, we set these hyperparameters so as to obtain optimal results 
in the sense of a trade-off between noise and actual discontinuities. 

The simulation image (see Fig. 3) is issued from insonification of the synthetic 
body shown in Fig. 1 (simulated transducer on top of the image). Attenuation 
is taken into account, and the pulse frequency is 7 MHz, which approaches, ac- 
cording to previous remarks, ideal conditions of imaging. Acoustic parameters 
of the tissues have been intentionally chosen so that the contrast of the envelope 
image is very poor. As a consequence, the result based on the envelope is unsat- 
isfactory: merely half of the initial contour is recovered. A different choice of the 
hyperparameters does not reveal new contours, but only artifacts, mainly due to 
attenuation. It is the introduction of the spectral contents through the complex 
reflection coefficient 71 that enables correct detection of actual boundaries, with 
relatively few artifacts. Therefore, in this particular situation, the use of the 
spectral contents appears to be necessary for correct boundary detection. 

The image of the human heart, acquired at 2 MHz over 106 acquisition lines, 
shows a long axis view of the left ventricle (see Fig. 4). Low acquisition fre- 




Fig. 3. Simulation image: (a) envelope, (b) boundaries detected based on envelope, (c) 
boundaries detected based on envelope with sub-optimal choice of hyper parameters, 
(d) boundaries detected based on 71 
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Fig. 4. Human heart: (a) envelope, (b) boundaries detected based on envelope, (c) and 
(d) boundaries detected based on 71 with two different choices of hyperparameters 



quency yields a poor visual quality of this image. The envelope enables to detect 
the internal walls of the cardiac cavities, as well as several cardiac structures. 
However, detection based on 71 gives completely random contours. Conclusions 
from section 4 are confirmed: at low acquisition frequencies, the spectral contents 
is irrelevant and there is little hope that it can be used to improve boundary 
detection. 

The image of the canine heart represents a long axis view of the left ventricle 
acquired at the frequency of 5 MHz over 120 lines (see Fig. 5). Higher acquisition 
frequency and large acquisition area result in a better image quality. It is worth 
noting that the structure below the left ventricle is a ghost (mirror) image, and 
therefore has no medical meaning. When using the envelope, an important part 
of the heart cavity is properly detected, along with different cardiac structures. 
Reflection coefficient 71 yields some new boundaries, such as closing the left 
ventricle at the apex. Reflection coefficient 72 yields similar results. In this sit- 
uation, the use of the spectral contents could contribute to accurate detection 
of contours. Nonetheless, parasite discontinuities appear, especially inside of the 
cardiac muscle where the spectrum shows very high variability. 




Fig. 5. Canine heart: (a) envelope, (b) boundaries detected based on envelope, (c) 
boundaries detected based on 71, (d) boundaries detected based on 72 
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6 Conclusion and Perspectives 

We have proposed to introduce the spectral contents of the RF signal to de- 
tect boundaries in echocardiographic images, and we have also described a novel 
method to perform this task. Results obtained on models, simulations and actual 
RF images converge towards the same conclusions. For low acquisition frequen- 
cies, the spectral parameters estimated from different tissues are not dissimilar 
enough to become meaningful for boundary detection. Indeed, the estimated 
parameters are subject to natural variability, due to fluctuations of the acous- 
tic properties inside one tissue, acquisition noise and estimation variability. In 
the case of higher acquisition frequencies, spectral parameters can successfully 
be used for detection of discontinuities, and we have shown a simulation case in 
which the introduction of the spectral contents significantly improves the quality 
of detected boundaries. We have also proven that complex reflection coefficients, 
in particular the coefficient of order 1, are suitable spectral parameters for the 
task of boundary detection. These conclusions open promising ways to improve 
boundary detection on echographic images. It appears feasible to use the spectral 
contents to find new boundaries, which are not detected based on the envelope 
image. It may also be possible to merge boundaries detected based on the enve- 
lope and several spectral parameters, providing a measure of confidence affected 
to the detected boundaries. 

The limitation on high frequencies is related to attenuation, which increases 
with frequency. Therefore, high acquisition frequencies cannot be used for imag- 
ing of the human heart, where distances involved are relatively long. Spectrum- 
based segmentation could be adapted for children (pulse frequencies up to 8 
MHz), and also transesophagus (transthoracic) echocardiography, where dis- 
tances are shorter. Besides, applications other than cardiac can be considered, 
such as intravascular imaging (pulse frequencies up to 40 MHz). 
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Abstract. The aim of the present study is to check, by means of ele- 
mentary mathematical tools, a conjecture according to which myocardial 
fibres are geodesic curves running on some surfaces. This conjecture was 
first stated and experimentally checked by Streeter (1979) for the equa- 
torial part of the left ventricle free wall. Quantitative polarized light 
microscopy provides measurements on fibre orientation that could lead 
to evidence that the conjecture remains true for the whole of the left 
ventricle. Study of the right ventricle is under progress. 

1 Introduction 

It is commonly believed [3] that the myocardium design and structure allow max- 
imal mechanical efficiency in the systole and diastole processes. The long-term 
purpose of our multi-disciplinary approach is to try and propose a model for the 
mechanical behaviour of the myocardium. In the long run, performance of the 
complete electro-mechanical system could be analyzed. For related works, see 
[1], [7]. It is well known that usual mechanical models for skeletal muscles are of 
no help for the myocardium. Obviously, the myocardium is not, as ordinary mus- 
cles, linked at both ends to a bone. Fibre micro-structure and fibre geometrical 
organization are quite different as well. We believe that the specific fibre orga- 
nization in the myocardium should be taken into account in a complete model. 
Numerous anatomical studies, see e.g. [8], [9], have been devoted along the years 
to a description of the fibres arrangement and we refer to [5] for an extensive 
bibliography. What will be sufficient to recall here is that the dissection or peel- 
ing techniques are not precise enough, since apparently preferred fibre directions 
can be inferred by the experimental process. Data have been improved by means 
of several techniques in microscopy, such as photonic or electronic microscopy. 
In the present work, we use the data provided by the polarized light microscopy 
devices developed by some of the authors [4] , [5] . 

More precisely, we intend to check the geometrical description proposed by 
Streeter [9] who introduced a topological representation of the left ventricle 
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as a “nested set of toroidal bodies of revolution” on which myocardial fibres 
run as geodesics. This information would increase our understanding of the 
biomechanical properties and propagation of electrical stimuli in the heart. 

The experimental technique, which is valid for both ventricles and for the 
septum, measures for a set of points located on several myocardial sections two 
angles from which the fibre orientation can be deduced. In other terms, the 
output of the experimental work is a discrete three-dimensional vector field. 
Assuming that the left ventricle has a structure of revolution, we use the invari- 
ance of Clairaut’s constant along geodesics as a first hint that the conjecture 
might be true. 

Work under progress is devoted to a similar description of the right ventricle. 
Note however that the simple structure of revolution of the left ventricle is no 
longer true for the right ventricle. 




Fig. 1. Maps of the azimuth and elevation angles obtained by means of polarized light 
microscopy in a coronal section. 



2 Data 

The data are obtained on fetal hearts. We refer to [4] and [5] for a complete 
description of the protocol. Let us just recall that the ventricles are embedded 
in a transparent resin in which they can be clearly seen after polymerization. 
Sections that can be transversal, sagittal or coronal can then be cut. A section 
thickness is 500/im, and, because of the thickness of the saw, adjacent sections 
are separated by a 250/im gap. The measurement technique relies on the birefrin- 
gence properties of myocardial cells: in short, the velocity of the light is slower 
when travelling along the long axis than along the short axis of the fibre. Results 
are given pointwise: each section is dicretized in 130/im x 130/im squares, and 
a mean angular information is collected for each of these elementary squares. 
For each voxel, the acquisition and representation processes result in two angles: 
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the angle of elevation 7e/e which is the angle between the fibre and the plane of 
the section, and the azimuth angle ^azi which is the angle between the projec- 
tion of the fibre on the section plane and a fixed direction in this plane, namely 
the east-west axis. A local knowledge at point (x, z) of these two angles which 
should range from 0 to tt for ^azi and — tt/2 to tt/2 for 7e/e determines completely 
the direction r of the fibre 



r = 



'^x cos 7eZe cos 7a2;i 

Ty = COS7eZeSin7a^i 

Tz = sin7eZe 



A drawback of the experimental device in its present state is that it provides the 
same value for 7e/e and — 7eZe- As a consequence and because of the averaging 
technique, values closed to 0° cannot be reached. Moreover, angles between 75° 
and 90° cannot be resolved. The accuracy of the method has been checked on 
myocardial samples which fibres are parallel: the resolution of the measurement 
method for both angles is 1°. 

From a mathematical point of view, we are given a discrete sample of a 
distributed vector field. 



3 Geodesics 

According to Streeter [9], in the equatorial part of the left ventricle free wall, 
fibres are organized into surfaces on which they run as geodesics. In order to 
check and extend this conjecture to the whole left ventricle, let us first recall 
some elementary properties of geodesics. 

Definition: A regular curve of a given surface of is a geodesic if and 
only if its principal normal at any point is normal to the surface. 

This definition has well known consequences in terms of lengths [2] . 

Proposition: Any geodesic locally minimizes the arc length between two 
points. Conversely, when a regular curve is the shortest path between any two 
of its points, it is a geodesic. 

Note that, in general, an arbitrary geodesic does not minimize globally 
the arc length between two of its points. For instance, an helix on a vertical 
cylinder links points on a same vertical line that can be more economically 
connected through a vertical straight line. Can we expect to find geodesics on 
any surface? Two answers are well known. First, local existence of a geodesic 
starting from a given point and tangent to a given vector belonging to the 
tangent plane at this point is an easy consequence of the fact that a geodesic 
on a parametrized surface is given by a system of two ordinary differential 
equations. Second, a global existence result is available in the case of closed 
surfaces: there exists a minimum geodesic joining two given points of a closed 
surface (particular case of the Hopf-Rinow theorem, [6]). 
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In the present study, we assume that the left ventricle has a structure 
of revolution. Namely, we see the outer surface as generated by rotating a 
meridian curve, described by Streeter as crescent-shaped, around a vertical 
axis, see Figure 2. The distance between this generating curve and the axis of 
revolution is close to 0, but not 0, at the apex which allows fibres to invaginate. 
Internal layers, if they exist, are also seen as rotation invariant. Computer 
Aided Design (CAD) approximate models using B-spline smooth surfaces have 
been constructed. They are out of the scope of the present paper and are 
used in a parallel work as test models for fibre reconstruction algorithms. This 
simplifying assumption allows us to use a specific property of geodesics of 
surfaces of revolution [2]. 




Fig. 2. An ideal setting. Left: nested meridian curves. Right: periodic geodesics deduced 
one from the other by rotations. 



Clairaut’s constant: For any point on a curve on a surface of revolution, let r 
be the distance from this point to the axis of revolution and let 0 be the angle 
between the tangent to this curve and the parallel of the surface at this point. 
If the curve is a geodesic, then the quantity rcosO remains invariant along the 
curve. It is called the Clairaut constant. 

Note that the above property is not sufficient for a curve to be a geodesic. 
Along all parallels, e.^., the quantity rcosO remains constant since it is equal 
to r, but only some of the parallels are geodesics. 

4 Experimental Validation 

As seems to be true from anatomical observations, assume that the left ventricle 
fibres are organized in nested toroidal layers with the same axis of revolution 
and that they are periodic geodesics of these layers, z.e., the fibres follow closed 
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curves with continuous tangents. Furthermore, assume that the fibres on a 
given layer can be deduced one from the other by rotations around the axis of 
revolution. Figure 2 describes an ideal setting where the geodesics have been 
computed on the abstract CAD model (the algorithm numerically solves the 
system of ordinary differential equations characterizing a geodesic). Then the 




Fig. 3. Left: isovalues of rcosO on a given coronal section. Right: Fibre trajectories 
crossing traces of a given C on several coronal sections. 



Clairaut constant is the same for all fibres in the layer, but may vary from one 
layer to another. As a consequence, if the conjecture is true, then the isovalues 
of r cos 6>, where 0 denotes here the angle between the vector r at point {x^y^z) 
and the parallel at this same point, will be nested toroidal surfaces. The traces 
of these surfaces on each coronal section (a section which is orthogonal to 
the left ventricle axis) should be concentric circles. This comes from the fact 
that the fibres are supposed to be organized in layers and globally rotationally 
invariant on each layer. Moreover, a same constant value of r cos 0 should appear 
on two separate concentric circles translating the fact that the fibre invaginates 
and, then, makes at least one complete turn before closing back. As for the 
trace of each surface of isovalues on a meridian plane, it should consist of two 
symmetrical closed crescent-shaped curves. 

With this in mind, we have computed r cos 0 for our experimental data. In 
order to use what is actually provided to us, namely the vector field r, we have 
used the obvious identity 

r cos 0 = \ — yTx + XTy \ . 

We have traced the lines of isovalues on several coronal section planes. We readily 
observe that on each section they are, as anticipated, concentric circles, and that 
a same constant value of r cos 0 appears on two separate concentric circles. In the 
same time, we have developed an algorithm which follows the fibres from section 
to section in the whole of the myocardium. By means of an interpolation pro- 
cedure, the algorithm first transforms the given discrete data into a distributed 
vector field. Then, for a given h, the iterative procedure constructs a point from 
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the previous one by moving with length h along its vector direction. In the left 
ventricle, the obtained trajectories cross the circles that correspond to a same 
Clairaut constant. In parallel, we have traced several isovalues of rcosO on two 
longitudinal section planes, one orthogonal and the other parallel to the inter- 
ventricular septum (namely, transversal and sagittal sections) , from the apex to 
the base. The result shows, as expected, a nested set of meridian curves which 
is similar to the abstract model described in Figure 2. 

Clearly, an arbitrary vector field r would not satisfay the above properties. 




Fig. 4. Traces of some values of C on longitudinal section planes parallel (left) and 
orthogonal (right) to the interventricular septum. 



The conclusion of this study is that checking the invariance of the Clairaut 
number against experimental values backs up the geodesic conjecture, at least 
for the left ventricle. In a future work, we intend to study the conjecture for the 
right ventricle. We will use the fibres tracking algorithm that has already been 
developed and used for the left ventricle. Simultaneously, an explanation of the 
geodesic structure in terms of mechanical efficiency will be investigated. 
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Abstract. Using a three-dimensional computer model of the human 
ventricular myocardium, we studied the role of ventricular architecture 
and conduction system in generating intramural activation patterns and 
the extracardiac electric field. The model represents the myocardium as 
an anisotropic bidomain; it incorporates detailed anatomical features, in- 
cluding intramural fiber rotation, the differences in the fiber arrangement 
of the trabeculae and papillary muscles, and a conduction system. Ec- 
topic activation was elicited at various depths, and “normal” activation 
was initiated via the conduction system. Extracardiac potentials were 
calculated throughout each activation sequence. The simulated epicar- 
dial potential maps resembled those measured in canine hearts, featuring 
a central minimum accompanied by two maxima in the early stages of 
ectopic activation, with the axis joining these extrema approximately 
parallel to the fibers near the pacing site. The simulated isochrones for 
the “normal” activation had characteristics very similar to those observed 
in isolated perfused human hearts. 



In the human ventricular myocardium, the anatomical arrangement of intra- 
mural fibers [19] and the gross network of interconnecting trabeculae on the 
endocardial surface, on which are numerous Purkinje strands, create a com- 
plex anatomical substrate for the propagation of electrical activity [1]. Both the 
activation sequence and, perforce, the intracardiac electrograms and/or body- 
surface electrocardiograms that are associated with it are determined by these 
anatomical features. 

1 Anisotropic Myocardial Structure 

The different physical properties observed along the axial and transverse axes 
of myocardial fibers (anisotropy) are a result of the geometrical arrangement 
of elongated cells and of the distribution of their gap junctions [21]. When an 
activation wave front initiated by point stimulation propagates through cardiac 
muscle, it assumes an ellipsoidal shape that reflects faster conduction along the 
fibers [5] ; the distribution of the potentials generated by the spreading wave front 
is also affected by the direction of the fibers. Experimental results obtained in 
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superfused myocardial laminas [22], on the epicardial and endocardial surfaces 
of isolated hearts [2] , and in a volume conductor surrounding an isolated heart 
[4] have shown that, in the initial stages of propagation, when activation elicited 
by point stimulation spreads through fibers that are approximately parallel, 
potential distributions show a central negative area surrounding the pacing site, 
accompanied by two positive areas, each with a maximum. The axis joining the 
two maxima is almost parallel to the fiber direction near the pacing site. 

The counterclockwise (CCW) rotation of fibers from epicardium to endo- 
cardium [19] further affects wave-front shape and associated potential distribu- 
tions in the later stages of activation. Taccardi et al. [23] recorded, with high 
spatial resolution, epicardial potential maps from the ventricular surface of ex- 
posed canine hearts during ventricular pacing. They identified features of the 
maps that reveal the direction of the myocardial fibers through which the ac- 
tivation is spreading and investigated how these features vary as a function of 
pacing site, intramural pacing depth, and post-pacing time. They observed that, 
during the later phases of propagated activation, the positive areas expand and 
rotate CCW when the activation wave front propagates from the epicardium to 
the endocardium, or clockwise (CW) when it propagates from the endocardium 
to the epicardium; the asymmetry of the potential distributions and multiple 
maxima appearing in the expanding positive areas were noted. 

Most of these results had been predicted by simulations of Colli Franzone 
et al. [3] and were substantiated by those of Henriquez et al. [9]. Both of these 
studies were based on bidomain theory and used a parallelepipedal slab of ven- 
tricular tissue with rotational anisotropy. Their simulated activation sequences 
initiated from epicardial, intramural, and subendocardial sites provided valu- 
able insight into the epicardial and intramural distributions of potentials and 
activation isochrones, but they were unable to replicate the fragmentation and 
asymmetry of the expanding and rotating positive areas. Simulating these char- 
acteristics would appear to require more realistic models. 

Many realistic models have been constructed to simulate propagated activa- 
tion in the whole heart [8]; among recent models, those introduced by Panfilov 
[20] and Berenfeld and Jalife [1] are of interest because they incorporate a rotat- 
ing anisotropy based on meticulously collected canine anatomical data [12,19]. In 
our previous study [11], we described the effects of the anisotropic structure on 
epicardial potential distributions in an anatomically realistic model of the human 
ventricular myocardium — based on previously developed propagation algorithm 
and calculation of extracardiac fields [14,18] — and we tested whether the features 
that could not be reproduced in idealized models [3,9] can be attributed to the 
anatomical characteristics of realistic structure. We compared simulated epicar- 
dial potentials for paced activation sequences with those recorded [23] during 
pacing at varying depths in the walls of both ventricles. 

Figure 1 shows epicardial maps for activation initiated at pacing sites within 
the right-ventricular (RV) free wall (5 mm thick, with fibers rotating from —44° 
at the epicardium to 32° at the endocardium). In these distributions, the major 
axis of the trough in the negative potentials was approximately parallel to the 
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Fig. 1. Simulated epicardial potential maps on a 17 x 33 grid at 10 ms after the onset 
of activation for pacing at different intramural depths in the right- ventricular free wall. 
The pacing sites were 0.5 mm apart, progressing from the epicardium (panel A) to 
the endocardium (panel J); the epicardial projection of each site is indicated by the 
filled circle. Isopotential lines are plotted for equal intervals, with no zero line; solid 
contours represent the positive and broken contours the negative values of potential; 
the magnitudes of the minimum and maximum are given (in mV) below each map. 
Note that the axis joining the maxima rotates counterclockwise with increasing pacing 
depths. From Hren et al. [11], with permission. 
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local fiber direction at the pacing site, and the positive regions rotated CCW 
as the pacing depth increased, following the transmural CCW rotation of fibers 
from the epicardium to the endocardium (Fig. IB-I). The distance between the 
maxima did not increase with the pacing depth, unlike in the thicker wall of the 
left ventricle. The maxima were unequal in strength, with the larger one often 
on the left side (Fig. lE-J) of the pacing site. 

Depending on the pacing depth, subendocardial and endocardial pacing ini- 
tially generated epicardial potential distributions that had a single oblong pos- 
itive area, which developed at 2-18 ms after the onset of activation into the 
usual pattern with one minimum and two maxima. Distributions calculated at 
10 ms after the onset of activation due to an endocardial stimulus at those pacing 
sites in the right ventricle where the free wall was thinner usually had a central 
negative area surrounded by a peripheral positive area with two maxima whose 
joining axis was approximately parallel to the local fiber direction (Fig. IJ). 
At subendocardial sites in a thick segment of the left-ventricular (LV) free wall 
with no trabeculae, on the other hand, as the pacing site approached the en- 
docardium, the distance between the maxima decreased, until, for endocardial 
pacing, a central oblong positive area emerged, from the major axis of which 
one could clearly deduce the local fiber direction. When the wall also had thick 
trabeculae, the positive area had no features that one could correlate with the 
local direction of the fibers. Interestingly, the magnitudes of the maxima did not 
change monotonically as a function of the pacing depth: they decreased with 
increased pacing depth until the pacing site approached the endocardium, then 
started increasing again (Fig. 1). The magnitude of the minimum, on the other 
hand, usually decreased monotonically with the pacing depth. This agreed with 
the results from the experimental study [23]. 

In general, our most important finding is that certain features of measured 
epicardial potential distributions — such as the fragmentation of the extrema and 
their asymmetric layout — which could not be reproduced by the slab models [3, 
9] , can be generated by a realistic model which features a variable wall thickness 
and recognizes that there are distinct differences in the fiber architecture of 
compact and trabecular regions of the ventricular wall. 



2 Specialized Conduction System 

Knowledge of the anatomy and distribution of the human atrioventricular (A- 
V) conduction system is important in understanding the sequence of “normal” 
ventricular activation and in establishing the anatomical basis of conduction dis- 
orders. The A-V conduction system begins at the A-V node and the bundle of 
His, proceeds through an area of fibrous tissue on top of the interventricular 
septum, and divides into the left and right bundle branches [6,17], which ram- 
ify subendocardially throughout the left and right ventricles into the Purkinje 
plexus. Fibers that fit the anatomical description of Purkinje fibers and show 
action potentials similar to those seen in the bundle of His are found throughout 
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Fig. 2. The geometry of the modelled conduction system in different views. The views 
of the LV septum (a) and LV free wall (b) show the complex structure stemming from 
the left bundle branch; note the vertical spiral gap (dashed area) in the LV conduction 
system in (b). The right bundle branch is shown on the RV septum (c) and the RV 
free wall (d). The Purkinje-myocardial junctions, at which the conduction system 
connects with the ventricular myocardium, are shown as light circles. 



the bundle branches, in false tendons [16], in an endocardial mesh, and in the 
chordae tendinae [10]. 

The anatomy of the left bundle branch (LBB) is extremely variable [25]. The 
LBB has usually been found to be initially a sheet-like structure, which then 
splits into three separate divisions with multiple interconnections [6] , then it fans 
out over the LV septal surface [17] and its divisions extend towards the bases of 
anterior and posterior papillary muscles, while the posterior rim is oriented to- 
wards the LV posterior papillary muscle. The right bundle branch (RBB) usually 
begins from the most distal part of the bundle of His [25], crosses over to the RV 
through fibrous tissue on top of the muscular interventricular septum, courses 
just beneath the endocardial surface of the RV septum before it emerges near 
the base of the RV anterior papillary muscle; from there it branches into fascicles 
leading to the posterior papillary muscle and the RV free wall. In the human 
heart, there is no evidence of penetration of the Purkinje fibers intramurally, 
as has been observed in other species [26]. From the Purkinje plexus, activa- 
tion is conveyed to the myocardium by transitional cells in Purkinje-myocardial 
junctions (PMJs) [1,26]. 

We have developed software for the interactive editing of the conduction sys- 
tem in our model of the ventricular myocardium. The user can create nodes 
on the endocardial surface by pointing to the desired locations, and connect 
them as required. These nodes can later be repositioned, and the connections 
can be either retained or modified. Some nodes are just bifurcation points for 
the conduction system; some are PMJs. The thickness of the connections and 
the extent of the PMJ patches are adjustable, to ensure connectivity and proper 
propagation between the conduction system and the myocardium. The propa- 
gation velocity in the conduction system was set to approximately 2.0 m/s. 
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Fig. 3. Simulated isochrones of activation in a three-dimensional computer model of 
the human ventricular myocardium. The horizonal layers — numbered 166, . . . , 26, from 
base to apex — correspond to sections depicted by Durrer et al. [7] ; each color represents 
a 5-ms interval, and the color coding of the activation times is the same as in [7] (except 
we use black to indicate the hrst 5 ms of activation in the conduction system). 
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Figure 2 shows one of our recent implementations of the conduction system. 
The initial network of the conduction system was guided by the anatomical 
descriptions in the literature [24,25]; subsequently, modifications were made in- 
teractively, to balance the timing of the initial activation. This was achieved by 
adjusting the length of the bundle branches, altering their course, varying the 
thickness of the segments between nodes, and altering the locations of the PMJs. 

The simulated activation sequence shown in Figure 3 agrees very well with 
isochrones constructed from experimental data obtained from isolated human 
hearts [7]. Very briefly, initial activity is on the left side of the septum (layer 
106), moving to the right, followed soon thereafter by inside-out activity in the 
apical and middle portions of both ventricles; the RV and LV breakthroughs 
take place at 25 ms and 30 ms, respectively. The areas that are activated last 
are in the basal RV near the pulmonary conus (layer 166). 

Modern techniques of clinical cardiac electrophysiology allow in situ investi- 
gation of the total excitation process in the human heart by recording electric 
signals from a large number of intracavitary [13] or endocardial [2] electrodes. 
These studies will provide an unprecedented impetus to the further development 
of computer models of the human heart. Conversely, models will be invaluable in 
synthesizing and interpreting the vast amount of information provided by these 
new investigative tools. This should lead to a better understanding of cardiac 
physiology and pathophysiology. 
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Abstract. In this study, the feasibility of two-dimensional strain rate 
estimation of the human heart in vivo is shown. To do this, ultrasonic 
B-mode data were captured at a high temporal resolution of 3.7 ms 
and processed off-line. The motion of the radio-frequency signal patterns 
within the two-dimensional sector image was tracked and used as the 
basis for strain rate estimation. Both axial and lateral motion and strain 
rate estimates showed a good agreement with the results obtained by 
more established, one-dimensional techniques. 



1 Introduction 

Regional strain and strain rate imaging have been introduced as new clinical 
tools to quantify regional myocardial function [1]. The most widely used method 
for cardiac strain rate estimation is based on myocardial velocity imaging, i.e. 
Doppler myocardial imaging [2]. This methodology has been validated both in 
vitro [3] and in vivo [4] and has been shown to be clinically applicable both 
in the experimental and clinical settings. However, a major drawback of the 
existing approaches is that they are limited to making only a one-dimensional 
measurement. Indeed, only the strain (rate) along an image line can be assessed. 
This limitation causes the technique to be angle dependent [5]. As a result, 
current ultrasound-based cardiac deformation data sets are incomplete. This 
could limit their application in quantifying cardiac function. 

In the field of ultrasound elastography other approaches to strain estimation 
by ultrasound have been described [6,7]. Some of these were shown to allow two- 
dimensional strain estimation by the tracking of radio frequency (RF) patterns 
within a two-dimensional RF image [8] . The purpose of this study was to test the 
applicability of such techniques to the heart and to make a multi-dimensional 
measurement of myocardial strain (rate) in vivo. 
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2 Methods 

2.1 Data Acquisition 

Cardiac ultrasound data were acquired from a young, healthy volunteer using an 
ultrasound scanner (GE Vingmed, System V, Horten, Norway) which allowed 
the continuous acquisition of digital IQ data (In-phase Quadrature sampled RF 
data [9]). A dedicated pulse sequence was developed in order to obtain B-mode 
data sets with high temporal resolution. To do this, the scan angle was reduced 
to 10 degrees and the number of RF image lines was set at 12. This resulted in 
an effective RF acquisition rate of 270 Hz. 

During apnea, a complete cardiac cycle of ultrasound data from the interven- 
tricular septum was acquired in an apical four chamber view using this dedicated 
pulse sequence. A standard 2.5 MHz phased array transducer (FPA 2.5MHz 1C) 
was used for fundamental imaging. Finally, the data acquired were transfered to 
a workstation for further off-line processing. 

2.2 Strain Rate Estimation 

RF data were reconstructed from the IQ data set. After axial and lateral linear 
interpolation of the RF data set (by a factor 2 and 8 respectively, resulting in spa- 
tial sampling of approximately 20 and 150 /im in the axial and lateral directions 
at the maximal image depth), the motion of the RF patterns between two suc- 
cessive RF frames was estimated using a newly developed algorithm which was 
based on the methodology described by previous authors [8]. However, in order 
to determine the time lag between two signal windows, the Sum of Absolute Dif- 
ferences (SAD) function [10] rather than the cross-correlation function was used. 
Moreover, temporal stretching (range : 0.97-1.03) of the first signal window, as 
presented in [11], showed to be indispensable. Window lengths were chosen to 
be 64 and 192 (interpolated) samples (corresponding to 1.2 and 3.7 mm) for the 
pattern to be tracked and the search region respectively. A window overlap of 
50% was used and the lateral search range was set to 17 (interpolated) image 
lines. 

The axial and lateral strains between two subsequent images were calculated 
as the axial and lateral spatial gradients of the axial and lateral motion esti- 
mates respectively [8]. Prior to applying a gradient operator, the axial motion 
estimates were median filtered (five pixels) to remove outliers and the lateral 
motion estimates were low-pass filtered by convolving them with a rectangular 
window of eleven pixels length. The gradient operator consisted of a linear fit 
through eleven consecutive motion estimates. Finally, the extracted strain values 
were normalized by the frame rate in order to obtain the axial and lateral strain 
rates respectively. 

3 Results 

Figure 1 (a) shows a M-mode image and an individual profile at a depth of 7 cm 
of the axial strain rate estimates over one R-R interval. The cycle starts with 
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a short period of negative strain rate which is followed by a longer period of 
muscle shortening, which ends approximately 260 ms after the onset of cardiac 
contraction. During the next 80 ms, strain rate values of different sign are found 
for the apical and basal regions. Subsequently, a bi-phasic lengthening is mea- 
sured, separated by a period in which length remains constant. The axial strain 
rates are homogeneous throughout the cardiac wall. 




Fig. 1. M-mode image and an individual profile of the axial (a) and lateral (b) strain 
rate estimates of the interventricular septum. The profiles were extracted at a depth 
of 7 cm. 



The same images are presented for the lateral strain rate estimates in fig- 
ure 1 (b). The lateral strain rate estimates are noisier than the axial ones but 
positive strain rates during the first 260 ms of the cardiac cycle can be observed. 
After this, a period of positive and negative lateral strain rates occur in the 
basal and apical regions respectively. Approximately 350 ms after the onset of 
the cardiac cycle, negative strain rates are measured throughout the cardiac 
wall. These are followed by a period of zero strain rate. Finally, another period 
of deformation is observed during the last 100 ms of the cardiac cycle. 

4 Discussion 

The strain rate patterns in normal individuals, derived by tissue Doppler meth- 
ods, have been described extensively [1,12,13]. The longitudinal strain rate of the 
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septal wall, as assessed from the apical four chamber view, is negative during 
systole (representing shortening of the muscle) and shows a bi-phasic lengthen- 
ing (positive strain rates) during diastole. The first diastolic lengthening phase 
is caused by a combination of relaxation of the cardiac muscle and early filling, 
while the second is a result of the increasing intraventricular pressure follow- 
ing filling due to atrial contraction. Both lengthening phases are separated by 
a period of constant length which corresponds to cardiac diastasis. It has been 
observed that isovolumetric contraction and relaxation are accompanied by a 
rapid change in shape of the ventricle [14]. This results in a significant amount 
of out-of-plane motion. Rapid changes in myocardial motion and strain rate es- 
timates have been reported during these periods. The average peak longitudinal 
strain rate values reported for systole, early and late diastolic filling are respec- 
tively : -1.5, 2.0 and 1.0 s“^ [13]. These have been shown to be homogeneous 
throughout the septum. 

To date, radial strain rates in the septum have not been studied by an ap- 
propriate method as current tissue Doppler methodology does not allow this. 
Indeed, radial strain rates can only be assessed from a parasternal transducer 
position. Using this approach, the septum is usually located in the near field of 
the probe. However, the radial strain rate of the posterior wall has been mea- 
sured accurately and it seems fair to assume that the septal wall shows a similar 
behavior : systolic thickening (positive strain rates) and bi-phasic diastolic thin- 
ning (negative strain rates). 

All of the regional strain rate characteristics already described in the liter- 
ature can be seen in figure 1 (a) and (b). In the apical four chamber view, the 
axial and lateral strain rate estimates correspond to the longitudinal and radial 
strain rates respectively [2]. Axial estimates were homogeneous throughout the 
septal wall; they showed shortening during systole and bi-phasic lengthening 
during diastole. Moreover, it can be observed that longitudinal deformation of 
the septal wall during atrial filling starts at the base and moves towards the 
apex. This phenomenon has well been described both in strain rate studies 
of the normal left ventricle and in computer simulations of left ventricular 
filling [15,16]. It is most likely the result of the propagation of a pressure wave 
and its associated blood vortex that enter the left ventricle through the mitral 
valve during atrial filling [15,16]. The radial strain rate estimates are noisier 
but the different phases of the cardiac cycle can still be identified : systolic 
thickening of the muscle followed by a bi-phasic diastolic phase. Both axial and 
lateral strain rate estimates had an isovolumetric relaxation time around 80 ms. 
This corresponds to the normal values in the literature [17]. 

Although our results demonstrate the feasibility of in vivo, two-dimensional 
myocardial strain rate estimation by means of RF tracking, much can be im- 
proved. Firstly, in order to construct the images presented in this report, no 
attempt was made to track the underlying tissue. This implies that the strain 
rate profiles shown, may contain information derived from differing anatomic re- 
gions. Nevertheless, the different phases of the cardiac cycle were distinguished. 
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Secondly, the lateral gradient operator was applied on all lateral motion esti- 
mates. As the outermost image lines were taken from the blood pool, this means 
that a possible source of noise in the lateral motion estimates was introduced. 
However, avoiding this effect requires more complex post-processing as the sep- 
tal wall moves laterally in the image during the cardiac cycle. Note that this 
lateral motion was clearly detected by the lateral motion estimator. This is one 
of the major advantages of this methodology : two-dimensional velocities, dis- 
placements, strain rates and strains can be extracted simultaneously. Finally, as 
the RF images were taken in a sector format, image lines diverge with depth. 
This implies that the intrinsic accuracy and noise characteristics of the lateral 
estimates will be depth dependent. Moreover, parameters such as line density, 
frame rate and interleaving will have a marked influence on the quality of the 
lateral motion and strain rate estimates. Tackling these problems in order to 
improve these two-dimensional estimates is left for future work. 

5 Conclusions 

In this paper the feasibility of two-dimensional strain rate estimation of the 
human heart in vivo was shown. Hereto, RF data sets of high temporal resolution 
were acquired and subsequently processed to extract motion and strain rate 
estimates based on a RF tracking algorithm. The properties of the extracted 
axial and lateral estimates showed a good correspondence to results found in the 
literature. Although many aspects could be improved, these initial results show 
that multi-dimensional myocardial strain rate estimation in vivo is feasible. 
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Abstract. This paper presents different ways to use the Doppler Tissue 
Imaging (DTI) in order to determine deformation of the cardiac wall. 
As an extra information added to the ultrasound images, the DTI gives 
the velocity in the direction of the sensor. We first show a way to track 
points along the cardiac wall in a M-Mode image (ID+t). This is based 
on energy minimization similar to a deformable grid. We then extend 
the ideas to finding the deformation field in a sequence of 2D images 
(2D+t). This is based on energy minimization including spatio-temporal 
regularization and a priori constraints. 

Keywords: Cardiac image processing. Ultrasound Image, Doppler Tis- 
sue Imaging, motion estimation, multi-modality image fusion. 



1 Introduction 

Doppler Tissue Imaging (DTI) is a recent technique which provides a partial 
information about the myocardial wall velocities. This is the velocity of the tis- 
sues in the direction of the sensor. This new data is represented by a colour 
added to the conventional ultrasound image. Two types of these images are used 
by the cardiologists. First 2D images as in figure I. a. are usually employed in 
ultrasound imaging. But the extra data which is collected during the acquisition 
implies a reduction of the temporal resolution compared to classical ultrasound. 
We cannot have more than a hundred images per second. M-Mode images (Fig- 
ure I.b) can give a better temporal resolution by reducing the studied area. They 
are obtained by choosing a ID segment on the image and watching its evolu- 
tion through time. Since there is only one dimension the temporal resolution is 
increased to about 500 frames per second. In this work, we used sequences of 
images of left ventricle of human hearts. 

Previously, M-mode images have been studied in [I] in order to track the car- 
diac wall with a variant of active contours. The contour C{t) = (t, /(t)) deforms 

to minimize E{f) = U Wl|/'(t)P + W2|/"(*)P + Pedge{C{t)) + Pvelocity(C{t))dt. 

Pedge ^ud Pyeiocity ^rc potentials which attract respectively C (t) to the edge and 
f'{t) close to the velocity measured by the DTI. [I] shows the advantage of using 
DTI in addition to the edge information. 

We will use different variational methods to study the deformation of the 
heart in M-mode and afterwards 2D DTI images. 



T. Katila et al. (Eds.): FIMH 2001, LNCS 2230, pp. 53-60, 2001. 
(c) Springer- Verlag Berlin Heidelberg 2001 




54 



V. Moreau, L.D. Cohen, and D. Pellerin 




(a) 2D Ultrasound Im- (b) M-Mode Ultrasound 

age. Image. 

Fig. 1. DTI Images. 




Fig. 2. Initialization of the steepest gradient descent. 



2 Tracking the Wall through M-Mode Images 



We now want to track several points through the cardiac wall on the M-Mode 
image making use of the velocity given by the DTI image 

Points are chosen for the initial time t = 0 along a hand given segment. They 
are regularly spaced on this vertical segment as seen on the left of figure 2. 

Let {Ci{t) = be the curves which perform the tracking 

of these points. If the points we want to follow are close enough, the curves 
{Cl,..., Cat} must be consistent with one another. In order to improve the 
tracking, we will consider the set of curves {Ci,...,Cat} as an active net: 
(^(t, s) = ys{t))s,t where only the second coordinate y{t, s) can vary. The active 
nets, or deformable grids, are defined in [2] and [3]. The net deforms according 
to the minimization of an energy. Our energy consists of three terms. 



1. A regularization term as in [2]. 

^regiy) = ff 0^(-^ ^ 



This term will be the one enabling interaction between successive curves. The 
first derivatives make the net contracts (which should be avoided by choosing 
a small) and the second derivatives enforce smoothness and rigidity. 








Deformation Field Estimation for the Cardiac Wall 



55 




(a) After a few iterations of steep- (b) Result of the tracking, 

est gradient descent. 

Fig. 3. Steepest gradient descent. 



2. An external term which attracts the derivative close to the given velocity 
measured by the DTI. This is an extension in two dimensions of [1]. 

Evelocityiy) = JUm - VDTi{t,y{t,s))'^dtds 

3. The last term is also an external term. We assume that the gray level is nearly 
constant along and around each curve. 

Egrayiy) = Jf s) + k) - l{0,y{0, s) + k)^dsdt 

We used both DTI and gray scale conventional ultrasound in this energy. 
The net y is obtained by the minimization of E{y) = Ereg{y) + i^Egray{y) + 
XEyeiocityiv) where A and y are postive constants, y must be chosen small be- 
cause Egray IS vcry sensitive to noise. We proceed as in [1] and [4] and use a 
steepest gradient descent, the discretization was done by finite differences. This 
gives in matrix form {Id + + tE(Y^) where Y is the discrete 

version of y and A is a sparse square matrix. To avoid a large matrix inversion, 
we can approximate {Id-\-rA)~^ ^ Id — rA for r small enough. We then apply a 
SOR algorithm [5]. The energy can have many local minima. Therefore we must 




(a) Optical Flow. Close up. (b) Velocity field obtained by 

Bottom of the lateral wall. our method on the same case. 

Fig. 4. Comparison between optical flow and the velocity field we obtained. In these 
images, the sensor is above and on the right of the pictures. 
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Fig. 5. Velocity field with the simplest method. 



use an initialization close to the solution. For this purpose, we choose an integra- 
tion of the velocity measured by the DTI: y{t + 1, s) « ^(t, s) ^(t, s))dt 

(figure 2) where dt is the time step. These curves, which are the direct interpreta- 
tion of the DTI velocity do not provide an accurate tracking of the wall. We can 
see that some curves leave the wall during the tracking. The hypothesis about 
the gray level is not exact in the case of ultrasound images. Nevertheless figure 3 
shows how this term of the energy tends to correct the tracking processed using 
only DTI. Figure 3.b presents a final result. We have compared these curves 
to manually traced curves by the cardiologist. We notice that the set of curves 
folow precisely the deformation of the cardiac wall. This automatic tracking has 
been of much help for cardiolgists in order to study various use of DTI images 

[6], [7]. 

3 Deformation Field in a 2D Image Sequence 

3.1 Presentation 

Our next interest is the study of 2D image sequences. The difficulties raised are 
the low temporal resolution and the incomplete DTI information. The DTI only 
measures the velocity in the direction of the sensor. Moreover the acquisition 
of the sequence instead of a single image implies a lower spatial resolution. Our 
motivation is to be able to diagnose a pathology with the help of DTI ultrasound 
images. For that purpose, we are mostly interested in recovering a complete 
deformation field. 
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3.2 Recovering the Velocity Field 

The DTI sequence contains two types of information: the velocity VDTi{x^y^t) 
in the direction of the sensor measured by the DTI and the conventional ultra- 
sound sequence I{x,y,t). We will use this second information to complete the 
first one. From the grayscale conventional ultrasound sequence, we can calculate 
the optical flow [8], [11], [9] and [10]. The optical flow is the distribution of ap- 
parent velocity of motion of brightness patterns in an image. The determination 
of the optical flow is based on the hypothesis that the brightness of a point 
I{x{t),y{t),t) is constant through time. Using the chain rule for differentiation 
we see that + ff = 0. If we let ^ and v = we have a 

linear equation in the two unknows u and v : IxU + lyV -\- It =0, where = |^, 




Fig. 6. Close up on figure 5. 




Fig. 7. Velocity field with a spatio-temporal regularization. 
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Fig. 8. Close up on figure 7. 




Fig. 9. Velocity field with a radiality constraint. 




Fig. 10. Close up on hgure 9 
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^ and It = This is called the constraint of the optical flow. It does 



not define a unique solution. Our solution is inspired by the Horn and Schunck’s 
method and will take advantage of the knowledge of the DTI information. The 
field we are looking for will satisfy three constraints: 

1. The optical flow constraint. 

2. The agreement with the DTI velocity. 

3. A regularity constraint. 

These properties are obtained as follows: we look for a vector field v) which 
minimizes 



E{u,v) = J j {xdtiu E Vdti'^ - '^DTi)‘^ E a{\\Vu\Y 



liv^lp) 



+ lyV + Itfdxdy 



where a and f3 are positive constants, Q is the image domain, {xdti^Vdti) 
denotes the direction of the sensor and vdti denotes the norm of the velocity 
measured by the DTI. We minimize this energy with a steepest gradient descent. 
We used a discretization by finite differences with an explicite scheme. We give 
an example in the figure 5. All the examples were obtained on the same image 
taken during the systole. The results are satisfying compared with simple optical 
flow (figure 4). We get a more regular field without increasing the diffusion. We 
observe that the velocity field mainly represents a motion of contraction. But 
this field is still very noisy. 

We can improve the results by using the idea of [11]. It consists of using a 
spatio-temporal regularization instead of a spatial regularization. It leads us to 
process the whole sequence simultaneously by minimizing the following energy: 



E{u,v) = [XDTIU + VDTIV - VDTif 

J J Jnx[o-,T] 

+f3{IxU + lyV + Ifj'^dxdydt. 



KllVt 



iiv^in 



The rest of the algorithm is the same. An example of result is shown if figure 7. 
As in [11], the result is much more coherent and complete. 

In order to improve results, we set soft constraints from a priori informa- 
tion about the heart deformation. In this type of image, the deformation is 
a contraction. Its center is the center of the heart. We then imposed a ra- 
diality constraint on the energy. Let O = {xq^Vo) be the center and ra- 
diuses Ri and R 2 define a ring around the myocardium wall: C{0, Ri, R 2 ) = 
{{x,y)^Ri < d{0,{x,y)) < R 2 }. We add a term to the energy of the form 
/[o-T]/ Ic(o Ri ~ ~ "^‘{y ~ yo))‘^dxdydt. We obtain figure 9. We can 

see the result of the constraint on the direction of the field but we could lose 
information while giving too much a priori information. The choice of adding 
this constraint or not is dependent on the application. 

We checked the validity of our results by deforming an image with the velocity 
field and comparing to the following image, both visually and with a quadratic 
norm. This comparison gave good results, but, since we calculate a smooth field. 
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it is not an quantitative evaluation of the quality of the velocity field. In future 
work, we will analyse the deformation field by the way of segmentation and 
classification algorithms. Then we will compare these analysis with other medical 
datas. That is how we will be able to check this deformation field estimation. 

4 Conclusion 

We have presented different ways to use the Doppler Tissue Imaging (DTI) in 
order to determine deformation of the cardiac wall. We first showed a way to 
track points along the cardiac wall in a M-Mode image (ID+t), based on a 
deformable grid energy minimization. We then showed how to estimate the de- 
formation field in a sequence of 2D images (2D+t), based on energy minimization 
including spatio-temporal regularization and a priori constraints. Future work 
includes analyzing the deformation of the cardiac wall and finding pathology 
starting from the deformation field method we described in this paper. 
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Abstract. Magnetic resonance tagging has proved to be an efficient 
non-invasive imaging technique for the study of the heart motion, pro- 
ducing a sharp and dense pattern on the tissues. Up to now, the quality 
of information derived from it has been unequalled by other modalities. 
It seems to be the perfect modality for the building of heart motion 
models. However, its main drawback was the tediousness of the image 
processing tasks it requires to extract the information. Recent achieve- 
ments [1,9] have overcome that. In previous work [2], we have presented 
a novel parametric class of deformation for the 2D heart motion in the 
short axis planes. From it, we can compute a very accurate and compact 
analytical expression of the walls motion from the grid information. We 
present here this kinetic model, and its applications to heart modeling. 
In particular, we build a decomposition basis of the 2D motion with only 
three orthogonal modes. 



1 Introduction 

1.1 Tagging Protocols 

Classical tagging, known as SPAMM (SPAtial Modulation of Magnetization) 
[14], produces a pattern (so-called tag)^ usually a grid, on the myocardial tissue, 
which remains along time, thus enabling to follow the motion. From Fisher’s 
improvement CSPAMM (for Complementary SPAMM) [6], which prevents the 
tags from fading away, Stuber et al [13] introduced a new method which en- 
ables to track the first acquisition slice along the cycle. This is known as the 
“Slice Following” technique. The true in-plane 2D motion of the initial slice is 
hence acquired, with few artifacts. Validations of this technique to healthy and 
pathological hearts have already been published [12,11]. 

* corresponding author 
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1.2 Motion Acquisition and Modeling 

In the last decade, many works have dedicated to the acquisition of relevant 
kinetic parameters from tagging [3,10,7,8,4,5]. All of these methods require to 
use 3D geometrical models to interpolate between the acquisition slices. This 
induces necessarily an additive bias, as the classical SPAMM images are highly 
anisotropic and subject to partial volume effects. The Slice-Following technique 
enables to avoid these drawbacks, and permits a direct “real” 2D motion recon- 
struction. 

In [1,2], we introduced a new algorithm which enables the fully automatic 
extraction of the tagging grid, in a very accurate and rapid manner. When 
applied to C SPAMM images, it even enables to extract the maximum intensity 
grid, multiplying hence the information by four. We then introduced a new 
parametric deformation class with which the grid motion in the short axis plane 
can be approximated with a very high sharpness. Moreover, the global motion 
parameters (rotation and contraction) are embedded in the model, and hence 
directly computed, for both walls. We briefly recall this deformation class, and 
then show an application of it: the approximation of the heart true 2D motion 
with only three orthogonal modes. 

2 Parametric Deformation Class 

2.1 Construction 

The motion of the normal LV walls in the short axis planes looks very close to a 
similarity deformation (rotation+contraction/dilation+translation). We chose to 
base our model on this transformation class. To make it closer to the reality, the 
idea is to enable the dilation and rotation to vary smoothly along the wall. Let us 
consider the complex plane centered at the LV barycenter (roughly computed). 
We introduce the deformation which transforms a point z = \z\e^^ into: 

N 

/(ai),d(^) = kl ^ , (a/c) G d G C ( 1 ) 

k=-N,k^0 

" V " 

^{ 9 ) 

\ai\ induces the homogeneous dilation/contraction, arg(ai) the rotation, d the 
translation, and all other Fourier terms stand for the “elastic perturbation” . The 
N order is the elasticity degree of the deformation in the “Fourier sense” . 



2.2 Computation of the Deformation from Tagging 

Let C be a contour (delineating a wall, the mid-myocardium, or any other line 
shaped zone of interest). We want to compute the deformation parameters of 
this contour, between times t — 1 and t. To do so, we weight each point of the 
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grid decreasingly with its distance to C at t — 1, and solve the Least Mean Square 
problem: 

= zt {weight{zt-i)) (2) 

With weight = T»aa:(i dist(z c)) ’ obtain an accurate interpolation scheme. As 
every point in the grid is taken into account, this scheme is also regularizing. 

2.3 Radial Strain 

A noteworthy property of our deformation class is that the radial strain is ex- 
pressed by: 

s^(0) = i^q0)i-i (3) 

where is defined in Eq. 1 for the deformation of contour C. 

2.4 Validations 

We have validated our interpolation class and scheme on data acquired from 15 
healthy volunteers, on a Philips Gyroscan 1.5 T [1,2]. The acquistion protocol 
was CSPAMM-Slice Following. The slice thickness was 8 mm, the grid spacing 4 
mm, and the temporal resolution 35 ms, resulting in sequences of 16 to 20 frames 
for a cycle. Acquisitions were taken at basal (10 mm below the valves), apical 
(10 mm above the radiological apex) and equatorial (equally spaced between 
apex and base) planes. 

We proceeded as follows: the endocardium and epicardium walls were manu- 
ally segmented on the first images, then the segmentations were propagated by 
applying the computed deformations. Visual validations of both segmentations 
and particular points following were considered extremely good. The parame- 
ter curves appeared highly accurate and better than those obtained by manual 
techniques [11]. The average absolute error on the fitting of the grid points was 
0.57 mm in radius and 1.06*^ in angle for the basal plane, and 0.45 mm and 1.12^ 
for the apical plane. This does not take into account the regularization effect of 
the deformation, over the error embedded in the grid. 

N = 3 appeared as the optimal order of the deformation, providing both 
sufficient elasticity and regularization. V > 3 is indeed “too elastic” for the 
noise level of the grid information, while N < 3 underfits the motion. 

2.5 Acquisition of Clinical and Kinematic Parameters 

As we dispose of an explicit analytical expression of the deformation at each 
phase, we can compute any parameter of clinical or modeling relevance (lumen 
surface, shortening profile, local motion amplitude, etc ... [12,11]). Moreover, 
the ai coefficient carries a huge kinetic sense: it provides, by construction, the 
best similarity approximation of the motion. Curves of temporal evolution of 
its modulus and argument are of great interest for both clinical and modeling 
purposes. See figures 1 and 2 for an example. 
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Fig. 1. Example of Tag Extraction Used for Validations (Basal) 



3 Towards a Compact Kinetic Model 



3.1 Complex Principal Component Analysis 



With A = 3, our deformation is guided by 7 complex parameters, which is rather 
compact related to the high degree of accuracy it handles. However, the main 
component of the motion is mainly controlled by the similarity coefficients ai 
and d, the rest being an elastic complement. We propose to reduce a posteriori 
this 5-parameter term to a single mode of variation on a cycle, for a same patient 
and same slice. Let Ct be the {ak)k=- 3 ...s vector at time t. Let us consider the 
classical hermitian product: < A\B >= A^B. If t = 0 is the reference time, 
Co = (0 . . . ai = 1 0 ... 0). Then < Co|Ct > Co = (0 . . . 0 ... 0) stands for 

the projection of Ct on the subspace of the similarities . Let now consider 



C,-<Co|C, >Co 
<Co|C, > 



( 4 ) 



Bt stands hence for the non-semi-rigid part of the deformation Ct, semi-rigidly 
registered to Co- 

A Principal Component Analysis (in C) of {Bt)t=o..tn leads to: 



Bt Mo+ < Vo\Bt > Vo 



( 5 ) 
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Fig. 2. Twist Computation for the Endo- and Epi-cardial Walls, at the Basal, Equa- 
torial and Apical Levels 




Fig. 3. Computation of Radial Contraction for the Endo- and Epi-cardial Walls, at the 
Basal, Equatorial and Apical Levels 
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Fig. 4. Left: Propagation of Segmentation (Basal, End- Systole). Right: Motion from 
Early to End-Diastole, at Basal Level, for Endocardium, Epicardium and Mid- 
Myocardium 



where Mq is the mean vector of and Vq the normalized first eigenvector. 
Calculations on our datasets of healthy patients proved that the mean vector 
Mo is always negligible, and that the first eigenvector has a very high weight 
related to the other ones. 

As, by construction, Bt is orthogonal to Co, so is Vb, and then we can rewrite, 
returning to vector C^: 



Ct -< Co|Ct > Co+ < K)|Ct > K) (6) 

Finally, by returning to Eq. 1, we obtain: 

3 

f(ai)Az] = aiz + d +A|^:| ^ , 9 = axg{z) (7) 

Similarity ^ ^ 

Elastic 

We then have an a posteriori orthogonal decomposition of the motion with 
three parameters: ai(t), d(t), A(t). 



3.2 Physiological Interpretation 

After preliminary results on our volunteers database, we have observed that the 
Vb eigenvector representing the elastic deformation had always similar shapes, 
its temporal influence being a unimodal curve, with maximum around the end- 
systole, and its spatial influence being maximal on three anatomical points: the 
two junctions of the right ventricle (RV), and the middle of the lateral wall. 
See Fig. 5. Physiologically, this can be interpreted by the fact that the septal 
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Fig. 5. Elastic Part of the Motion, for Both Walls, at End-Systole, Basal Shot. Scale 
*4 for Visibility 



wall has a diphasic motion, taking part in both the LV and the RV motions, 
while the lateral wall is almost free of influence from the RV. In a sense, the RV 
perturbates the LV motion from a perfect similarity. 



4 Conclusion 

We have presented a new parametric deformation class from which the 2D motion 
of the Left Ventricle in the short axis planes can be very sharply approximated 
and modeled. Coupled with the “CSPAMM-Slice Following” tagging technique, 
it enables to compute the true motion of the left ventricle in the short axis plane. 

By analyzing our deformation parameters, we have showed that the non-rigid 
part of the motion can be reduced to a single-mode deformation, which can be 
physiologically interpreted. 

The three-orthogonal-mode decomposition of the motion we obtain will be of 
great help in many image processing tasks, including 2D registration and motion 
computation in non-marked images, as also in model-guided segmentation. 

Our method was very recently extended to 4D modeling, and very promissing 
results have already been obtained. 
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Abstract. An image processing pipeline is presented for the quantitative analy- 
sis of 2D grid tagged Magnetic Resonance images. The first step concerns the 
automatic extraction of the tagging pattern and the definition of the left ven- 
tricular myocardial contours. In a second step, a spatio-temporal displacement 
field is fitted to the tag data points. Finally, parameters related to the contractile 
function can be investigated through graphic displays, movies and statistical 
analysis. 



1 Introduction 

Magnetic Resonance (MR) tagging is a technique of choice for analyzing the motion 
of the heart. Since the introduction of the technique more than ten years ago [1], a 
great number of methods have been published for the processing of 2D and 3D tagged 
images [2-10]. Up to now, despite the availability of MR tagging sequences on most 
MR scanners, no standard tool has been commonly adopted for the quantitative analy- 
sis of the heart motion from MR tagging. In fact, the quantitative analysis of tagged 
MR Images requires several successive steps (Figure 1) each of which having to be 
carefully validated. In a previous work, we have proposed a simple approach for fit- 
ting a spatio-temporal displacement field [10] from which several parameters related 
to the mechanics of the myocardium can be derived. Some tools have also been devel- 
oped for the analysis of parameter sets in order to detect contraction abnormalities 
[ 12 ]. 

In this paper, we describe the processing scheme, named TagAnalyze, for the com- 
plete analysis of 2D + time tagged data from the acquisition to the computer aided 
diagnosis of contractile abnormalities. 



2 Processing Pipeline of Tagged MR Images 

TagAnalyze is a processing pipeline for 2D tagged MR Image sequences. It handles 
grid tagged images in order to extract clinically relevant parameters. Figure 1 illus- 
trates the main steps of the processing pipeline, which are detailed in the subsequent 
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subsections. Three examinations were considered in this study : one healthy volunteer, 
two ischemic patients at rest and under pharmacological stress. 



Subject 

i 

Tagged MR Image Acquisition 



I 

Tagged MR Image Processing : 
Tags and contours extraction 
(semi-automatic step) 

Deformation Computation 

Deformation Analysis using FDA 

Computer Aided Diagnosis 
of the Contractile Function 




Fig. 1. Main steps of the proeessing pipeline of the tagged MR images from aequisition to 
diagnosis. A short axis sliee at six instants through the eardiae systole (a-f) is shown on the 
right. FDA stands for “funetional data analysis”. 



2.1 Tags and Contours Extraction 

Our approach requires the extraction of the tagging pattern and of the left ventricular 
(LV) contours. First, a Region of Interest (ROI) is interactively defined. A cine loop of 
the image sequence allows the user to check the appropriateness of the ROI location 
and size. Then, a myocardial mask is defined by the manual tracing of the epicardial 
contour onto the end-diastolic image and of the endocardial contour onto the end- 
systolic image. Extraction of the tags and displacement calculation will be performed 
within the region given by the mask. 

Extraction of the tagging pattern is based on the concept of active meshes [13-15]. 
The initial position of the grid tagging pattern is interactively defined. Each tag direc- 
tion is defined one after the other by tracing a line onto the image. Intersections 
(nodes) of the grid are automatically derived and a mesh structure is sent to the tag 
extraction algorithm (Figure 2). Each node is assimilated to a particle of mass m 
linked to its neighbors with springs (stiffness k). The spring net, initially at rest, is 
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submitted to an external force field originating from the tagged image. An equilibrium 
is reached when the following global energy is minimized : 

MU + CU + KU ^F{t) 

where M is the mass matrix, C the damping diagonal matrix (damping force is propor- 
tional to the velocity and defined by the constant c), K is the stiffness matrix (the stiff- 
ness of the springs is set to k), F is the force matrix and U the matrix of the node dis- 
placements. The set of linear equations is solved using finite differences and an itera- 
tive process. The external force field is computed from the image itself by a Canny- 
Deriche gradient operator. Note that nodes outside the myocardial mask are fixed by 
setting their mass to infinity while nodes that are inside have a unit mass. Also, the 
convergence process is interruptible so that interactive corrections are made possible 
(automatic process can be resumed after manual editing). The resulting mesh at the 
current frame is used as the initial solution for the next frame. 




Fig. 2. (a) Definition of the initial regular mesh, (h) Resulting mesh at eonvergenee on the first 
frame (end-diastole) for the healthy subjeet. 



The tags’ location is further refined using an active contour which is a simple and 
fast geometrical model similar to the discrete dynamic contour model of Lobregt and 
Viergever [16]. It is applied individually on each tag line which is first subdivided and 
constrained in the vicinity of the original segments (Figure 3). Usually, 30 iterations 
are sufficient and the process rapidly ends up when the speed of the contour points is 
close to zero. This step is important because each tag line is defined by a limited num- 
ber of points within the active mesh. The smoothing step adds points to the tag lines 
and regularizes their shape while constraining them to lie inside the tag valleys. The 
consistency through space and time of the tags is verified with a cine-loop. At any 
stage, user interaction is allowed. The tag points within the myocardial mask for each 
frame are finally saved in a file. 
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Fig. 3. Smoothed tag lines using the diserete dynamie eontour model in the first three frames of 
the healthy subjeet. 

2.2 Displacement Field Estimation 

A spatio-temporal displacement field is fitted to the tag points using the algorithm 
described in [10]. Briefly, the components of the inverse displacement field (from 
deformed state to the reference undeformed state) are first estimated by computing the 
parameters of a hierarchical cosine series model. Using the correspondences given by 
this model, the parameters of the direct displacement field are calculated. Optimal 
model orders have been derived from an accuracy study as a function of the tagging 
pattern parameters [10]. 

One advantage of such a parametric model is that it allows for the straightforward 
calculation of motion related parameters such as deformation tensors, principal defor- 
mations, speed and acceleration. The mean global translation is subtracted in order to 
keep only the deformation part of the myocardial media. Displacement and deforma- 
tion tensors are displayed using t\\Q fur-like representation based on 2D texture synthe- 
sis [17] where the fibers encode the orientation, direction, and magnitude (color) of 
the vectors (Figure 4). 




Fig. 4. Representation of the displaeement field at end-systole for patient P2 at rest with the fur- 
like representation (spatial orders = 5, temporal order = 7). 
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The concatenation of the representation at each time point of the cardiac cycle al- 
lows for the visualization of the temporal evolution of the parameters. 

2.3 Analysis of the Contractile Function 

The previous step provides a wealth of information that needs to be summed up for the 
clinician. Evolution of deformation parameters can be displayed using multi-box plots 
and bull’s eye representations. Figure 5 shows the temporal evolution of the maximum 
shortening at midwall and for the medium level with the healthy subject and the 
pathological case PI. Statistical tools have been developed in order to ease the analysis 
of the parameters through examinations. In particular, a healthy subject database has 
been studied to derive a statistical model of the normal contraction pattern from which 
a given patient can be compared to. Interested readers are referred to [1 1-12]. 
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Fig. 5. Evolution of the Maximum Shortening at a medium level for the healthy suhjeet (up) and 
patient PI at rest (bottom). 



3 Implementation Considerations 

The programs related to the tagged image processing and the displacement field fitting 
have been written in C-C++ language and embedded in a Tcl-Tk^ user interface. Some 
graphical functions make use of the Vtk^ library. It is therefore platform independent 
and do not require any commercial libraries. 

The automatic steps (active mesh and contours) of the tags and contours extraction 
from one 2D tagged sequence are very fast (a few seconds). User control is still re- 
quired especially at the mask borders. Less than 3 minutes are needed for the spatio- 
temporal model estimation with the optimal orders (Spatial orders = 5, temporal order 
= 1,9 phases) including the generation of the displacement field representations on a 
Pentium III, 800MHz. 



^ Ajuba Solutions 
^ Visualization ToolKit. 
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4 Discussion and Conclusion 

We have presented the 2D tagged image processing pipeline, TagAnalyze. Two di- 
mensional grid tagged images can be fully exploited within a reasonable amount of 
time. The extraction of the tagging pattern is fast but still requires user control. We 
have insisted on the importance of the tag smoothing step. Its impact onto the dis- 
placement field fitting could be evaluated. Note that multi-slice short axis acquisitions, 
covering the LV, can be treated by the pipeline as well. 

The remaining important problem of this pipeline concerns the automatic extraction 
of the endocardial and epicardial contours which was manually performed by the 
definition of a mask. The mask edges need not be accurate for the displacement field 
computation. In practice, we usually overestimate the LV epicardial border while 
underestimating the endocardial border. Although not crucial for the calculation of the 
myocardial deformation, accurate definition of the LV contours is required for the 
interpretation of the results. The automatic LV segmentation is a non trivial task since 
the tags blur the contours. If no standard cine acquisition is available, one needs to 
preprocess the images in order to clear up the tags. Denney overcomes the problem of 
myocardial contour identification introducing the ML/MAP method for the identifica- 
tion of tags without user defined contours [8]. It is based on the physics of the tagging 
process and on signal intensity and noise statistics computed from image data. The 
method has been applied to parallel tag patterns. In our method, only a crude estimate 
of the inner and outer LV contours is used to define a ROI. The MAP hypothesis test 
introduced in [8] could be used to improve a posteriori the definition of the contours at 
each phase. 

An original heart motion estimation method that uses harmonic phase (HARP) im- 
ages as been proposed by Osman et al. [9]. This method relies on basic mathematical 
operations and is therefore fast. It allows for the reconstruction of displacement fields 
for small motions. Comparison of both approaches would be particularly interesting. 

The extension of the pipeline to process 3D tagged images implies the adaptation of 
the pipeline for the treatment of long axis images. The checking of the consistency of 
the 3D tags in space and time would necessitate 3D visualization tools. Based on the 
same principle, the extension of the displacement field computation algorithm has 
been performed but, since the process has not yet been optimized, the 3D computation 
is time consuming. 

TagAnalyze is currently involved in the study of patients suffering of hypertrophic 
myocardial diseases. 
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Abstract. Methods for the precise measurement of three dimensional 
myocardial motion non-invasivley with MRI have recently been developed. 
These methods use a technique called “presaturation tagging” to mark the 
myocardium, and rapid MRI to track the motion of these markers. A unique 
capability of this method is the production of strain images representing the 
local deformation of the myocardium. These images clearly show the sequence 
of events during the activation of the heart, and can demonstrate abnormalities 
caused by asynchronous electrical activation or ischemia. Coupled with the near 
simultaneous mapping of electrical depolarization with an epicardial sock array, 
we can investigate the relationship between electical activity and mechanical 
function on a local level. Registered fiber angle maps can be obtained in the 
same heart with diffusion MRI to assist in the construction of the mechano- 
electrical model of the whole heart. 



1 Introduction 

The ability to measure the precise mechanical function, the electrical function, and the 
underlying fiber architecture in the same heart in vivo may uncover the interactions of 
these constituents in normal and abnormal cardiac function. 

The measurement of the electrical activity of the heart is a mature field of research 
employing intra-cavity electrodes [8], baskets, optical techniques [13], monophasic 
action potentials [5] and body surface electrode mapping [2,3,9] among other 
techniques. 
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Methods for measuring local three dimensional myocardial motion non-invasively 
with MRI have been developed using presaturation tagging patterns [1,10,17]. 
Recently, methods for measuring the diffusion tensor in vivo have lead to a method for 
measuring the fiber angle in soft tissue [6,14]. 

An experimental protocol has recently been developed in which electrical mapping, 
myocardial strain mapping and fiber angle mapping can be achieved in the same heart 
[4]. The data are registered so that local correlations can be made between these three 
features. 



2 Cardiac Tagging Techniques 

In cardiac tagging a set of saturation pulses placed in the tissue provides a signal 
intensity pattern in the tissue; the change in shape of the intensity pattern in the image 
reflects the change in shape of the underlying body containing the intensity pattern. 
Originally demonstrated by Zerhouni et. al. [17] with saturation pulses, and by Axel 
with SPAMM pulses [1], it was shown that parallel lines can be used to mark the 
tissue effectively. 

The objective of the analysis of tagged images is to track the 3D motion of each 
material point in the heart, and then to compute the six components of the strain tensor 
at each point for a sequence of time points throughout the heart cycle. The strain 
tensor characterizes the local deformation of the myocardium. Bulk translations and 
rotations of the entire heart may actually dominate the displacement and velocity 
measurements, but these are of limited value as an index of local myocardial 
contraction. In order to obtain precise quantification of the regional strains, the 
position of the tags must be measured with a “tag detection” algorithm [7]. Once the 
relative position of the tags have been determined as a function of time, these data can 
be used to estimate the strain tensor at each point in the myocardium. One method for 
doing this is a displacement field model based on B-splines [12]. 



3 Measurement of Myocardial Function During Asynchronous 
Activation 

In order to evaluate the relationship between electrical excitation and the onset of 
mechanical contraction, MR tagging experiments were performed during ectopic 
pacing in anesthetized normal dogs [4,11,15]. When systolic contraction was evoked 
by right atrial pacing, the LV was excited via the normal pathway of the Purkinje 
system and the pattern of mechanical activation was found to be very uniform as a 
function of position. However, when the heart was paced from a ventricular site, 
significant asynchronous and spatially heterogeneous contraction was observed. 

The precise sequence of events during ectopic excitation was particularly evident 
on the strain images. Fig. 1 shows the evolution in time of the circumferential 
component of the 3D strain tensor (Ecc) evaluated at the mid-wall for two pacing 
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sites; because the fiber direction is essentially circumferencial at the midwall, this 
component of the strain tensor closely matches muscle fiber shortening. During atrial 
pacing the ventricle is activated through the normal Purkinje pathway and muscle 
shortening evolves relatively homogeneously over the ventricle; this is demonstrated 
with the uniformly increasing blue color over the ventricle in the top row of Fig. 1. 
With epicardial ventricular pacing, early mechanical activation was observed at the 
pacing site, followed by propagation of a contraction wave-front to the opposite side 
of the heart. This is seen as the blue "wave" of muscle shortening emanating from the 
RV freewall pacing site in the second row of Fig. 1. There was also a significant 
"prestretch" of the late activated myocardium on the opposite side of the heart from 
the pacing site, shown as a bright yellow color. This prestretch was quite pronounced 
(15-20% in some cases) and occurred in the first 100ms after the ventricular pacing 
pulse. 




Fig. 1. These eolorized surfaees show the eireumfereneial streteh and eontraetion of the midwall 
of the myoeardium; the RV extends from 5 o’eloek to 10 o’eloek. This surfaee runs through the 
FV and RV. The top row shows the initial streteh during filling (yellow eolor) and a uniform 
eontraetion over the entire heart (blue eolor). The bottom row shows the same heart during 
eetopie paeing on the epieardial surfaee of the RV. The FV freewall udergoes a streteh during 
RV eontraetion, eontraeting after a delay of approximately 100ms. 

An alternative way of visualizing the contraction pattern is to graph the time course 
of strain for each material point of the LV. Each graph can be mapped to a position in 
an array that corresponds to a position in the LV. An example of such a LV strain 
map is shown in Fig. 2 where the sequence of circumferencial stretching and 
shortening (mid-wall Ecc) for selected midwall sites in the LV freewall and the RV 
freewall. The points shown ar subsampled from a larger array of data. The LV strain 
evolution maps are an excellent method for observing the delay to the onset of 
contraction for different regions of the heart. 

While asynchronous mechanical activation is obvious from both the color encoded 
strain images shown in Fig. 1 and the LV strain map showing Ecc vs. time in Fig. 2, 
we can also define the “mechanical activation time” as the time at which the muscle 
begins to shorten. This will correspond to the time at which the prestretch or 
ventricular filling reaches a peak and myocardial shortening begins, as shown by the 
red vertical lines in the LV strain maps of Fig. 2. 
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Fig. 2. Eaeh box represents the eireumfereneial strain vs. time over a 400 ms time period during 
prestreteh and eontraetion. The top row is the base of the heart, the bottom is the apieal region. 
The eolumns represent (from left to right) septum, anterior wall, lateral wall and inferior wall. 
The red vertieal bar represents the “meehanieal aetivation time” for eaeh loeation during RV 
paeing. 




The mechanical data described above can be obtained from a dog in vivo while 
nearly simultaneously measuring epicardial electrical excitation maps with an 
electrode array placed around the heart. (The electrical recordings are obtained 
between MRI image acquisitions.) From these data, maps of both mechanical 
activation (time of onset of shortening) and electrical activation (time of maximum 
negative slope in unipolar recordings) were computed. Voltage vs. time plots for the 
same locations as shown in the Ecc plots of Fig. 2 are shown in Fig 3. 



4 Fiber Angle Measurement and Modelling 

Fig. 4a and b show maps of fiber orientation for the heart measured in vitro with 
diffusion tensor imaging. These measurements are made after fixation in order to 
achieve very high resolution (0.7 x 0.7 x 0.9 mm voxels) using long imaging times 
(days). The diffusion of water in the muscle was sampled with diffusion encodings 
oriented in 28 different directions, each of which had a diffusion length of ~700nm. 
The map shows a clear transmural gradient in fiber angle from the epicardium to the 
endocardium. The raw fiber angle data can be used to feed a finite element model to 
describe a smooth field of fiber angles over the heart as shown in Fig. 4b. After the 
computation of the principle directions of diffusion each eigenvector is projected to 
the surface to compute fiber angles. These fiber angles are fitted to bilinear hermite 
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Fig. 3. Unipolar voltage measurements taken from the surface of the heart with a sock electrode. 
The locations are the same as those shown in Figure 2. The red asterisk represents the time 
chosen as “electrical activation time” during RV pacing. 




Fig. 4a. A “fiber angle map” in a short axis plane. The color represents the z component of the 
principal eigen-vector from the diffusion tensor at each location. Blue means the vector is in 
the plane, red means the vector is orthogonal to the plane. (Data courtesy of Pat Helm and Rai 
Winslow) 
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Fig. 4b. A finite element model showing the direetion of the myoeardial fibers as modeled from 
the diffusion tensor data. (Data eourtesy of Pat Helm and Rai Winslow) 

functions over the volume of the heart. The fine structure of this mapping will be key 
in determining the activation pathways for the propagation of electrical activation, and 
it may also have a bearing on the mechanical performance of the tissue. 



5 Discussion 

Strain imaging derived from MRI tagging data combined with electrical mapping 
using an electrode sock now gives us a new tool for studying the relationship of the 
temporal kinetics of the electrical and mechanical function during myocardial 
contraction. High resolution, co-registered fiber angle maps in the same heart will 
also allow us to investigate the underlying architectural substrate of this function. 

It is now possible to study the effect of mechanical stretch on the electrical nature 
of the heart in vivo. It has been demonstrated that rapid prestretch will stimulate 
depolarization of the myocardium especially in regions that have the greatest 
compliance and experience greater relative stretch [5], and that the timing of the 
application of prestretch determines if an arrhythmic depolarization of the LV will 
occur [16]. This new strain imaging technique combined with simultaneous electrical 
mapping will allow us to investigate the local prestretch needed to generate these 
arrhythmic beats. 

Acknowledgments. Many of the results and conclusions reported here are from the 
collaborative efforts of the Cardiac MRI Research Group at Johns Hopkins and the 
Medical Imaging Group at NHLBI. The authors would especially like to acknowledge 
the efforts of Scott Chesnick, Michael Guttman, Chris Moore, Walter O’Dell, Scott 
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Abstract. In this paper, an approach for the assessment of 3-D func- 
tional maps of the heart is proposed. It relies on the model-based co- 
registration of MR anatomical and PET metabolic images and the ex- 
traction of an individualized anatomical heart model from MR images. 
This results in a 3-D geometrical model of the heart for which functional 
parameters such as FDG uptake can be attributed and visualized. 



1 Introduction 

Cardiac viability studies mainly rely on the integration of the myocardial per- 
fusion, metabolism and contractile function. Metabolism can be analyzed with 
Positron Emission Tomography (PET) and contractile function can be stud- 
ied using Magnetic Resonance Imaging (MRI), for instance. In order to assess 
the cardiac viability assessment, it is crucial to accurately co-register the multi- 
modality images in order to be able corroborate the complementary functional 
information. In the “MunichHeart” software, endocardial and epicardial con- 
tours are manually delineated in Short Axis (SA) MR images and coregistered 
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with the same contours extracted from PET or SPECT (Single Photon Emission 
Tomography) using the maximum count detection algorithm [1]. In [2], maxi- 
mal myocardial deformation from tagged MRI and EDG-PET metabolism was 
combined using fuzzy rules to generate polar maps representing the viability. 
The registration was performed manually with the help of the long axis angles 
defined in the MRI protocol. In this paper, we present a model based approach 
for the automatic registration and segmentation of the right (RV) and left ven- 
tricles (LV) of the heart in PET and MR images. It provides a 3-D geometric 
representation on which functional information can be displayed. The data and 
the approach are presented in Subsections 2.1 and 2.2, respectively. Results are 
presented in Section 3 and discussed in Section 4. 

2 Material and Method 

2.1 Cardiac Imaging Protocol 

The data set is composed of MR and PET images of ten patients suffering from 
three vessels coronary artery disease [3]. Mean age was 69 (8 men, 2 women). 
All patients underwent MR and fluorine- 1 8- deoxy glucose (EDG) PET imaging 
within 10 days. 

The MR imaging was performed at the Department of Radiology of Helsinki 
University Gentral Hospital with a 1.5 T Siemens Magnetom Vision imager 
(Siemens, Erlangen, Germany). A series of 39 EGG-gated contiguous transaxial 
images was acquired during free respiration using TurboELASH sequence (Eig. 
la). The pixel size and the slice thickness were 1.95 x 1.95 mm and 10 mm, 
respectively. Eive EGG-gated breath- hold cine SA sections 15 mm apart were 
also acquired (Eig. lb). The pixel size for SA slices was 1.25 x 1.25 mm and 
slice thickness 7 mm. About 15 time points were taken for each section with the 
repetition time of 40 msec. 

Static PET imaging was performed using a Siemens EGAT (Siemens/GTI, 
Knoxville, USA) PET scanner. A series of 16 contiguous transmission and emis- 
sion images was acquired. The pixel size and the slice thickness were 2.41 x 2.41 
mm and 6.75 mm, respectively (Eig. Ic, d). 




Fig. 1. (a) Transaxial and (b) SA MR images, (c) transmission and (d) emission PET 
images. 
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2.2 Method Overview 

The aim of the overall approach is to extract a 3-D anatomical model of the 
heart from patient MR images and to incorporate functional data, such as FDG 
uptake and other clinically relevant parameters to the model. The process is sum- 
marized in Fig. 2. First, MR transaxial images are co-registered with the PET 
transmission image. The obtained registration parameters are used to register 
transaxial MR images and the PET emission image. PET images that corre- 
sponds SA MR images are calculated by using MR header information. Then, 
a 3-D biventricular model is initialized in the SA MR images to segment the 
myocardium. The deformed model is finally transformed into the registered SA 
PET image to obtain EDG uptake values. The main steps of the method are 
summarized in the following sections. 




Fig. 2. Extraction of the 3-D anatomical and functional model of the heart including 
PET-FDG uptake : main steps of the process. 



2.3 PET-MR Image Registration 

The method for registration of MR and PET images is fully described in [4] and 
[5] . It minimizes the distance between the point set from segmented PET trans- 
mission image surfaces and a distance map [6] built upon the segmented transax- 
ial MR surfaces. The thorax and lungs surfaces are segmented from transaxial 
MR and transmission PET images by using the deformable model based seg- 
mentation method presented in [7] . S A PET slices are computed from registered 
images by using MR header information (Eig. 3). 



2.4 Segmentation of the SA MR Images 

We use a 3-D prior biventricular deformable model of the myocardium to simul- 
taneously segment the LV and the RV in SA MR images. In order to constrain 
the segmentation, we introduce in our model a priori information about topol- 
ogy, geometry (Eig 4. a, b) and physics of the heart muscle, considering it as a 
linear elastic medium. 

Eirst, the template is manually initialized in the MR volume (Eig 5.). Then 
it is submitted to a force field derived from the image and elastically deformed 
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Fig. 3. The 2-D representation of registered (computed) end-diastolic SA MR (top) and 
PET emission (bottom) image slices for the HI case. In the middle row, the images are 
overlaid in block format to visualize registration results. 



until its boundaries fit the contours extracted from the data. We thus obtain a 
3-D geometrical representation of the patient’s heart. Fig. 6 shows segmentation 
results obtained for 2 patients. This approach has been fully described in [8]. 

2.5 Extracting Functional Data Associated to the Model 

Functional data can be attributed to the cells and nodes of the 3-D biventricular 
model. Some parameters with clinical interest can be derived from MR Imaging. 
First, the model is labeled so that it is possible to directly separate the right and 
left ventricles or obtain internal and external surfaces. Moreover, cavity volumes, 
myocardial mass and local wall thickness are computed. In order to enrich the 
model with metabolic information the model is transformed into the registered 
PET-FDG emission image. Medial surface is automatically calculated between 
LV endo- and epicardial surfaces. The medial surface nodes, M^, are calculated 
as follows. For each node Ai, of the endocardial surface, we compute the normal 
to the surface that intersects the epicardial surface at Bi. Mi is taken as the 
middle point of the line [AiBi] (Fig. 7a). Surfaces are transformed to registered 
PET emission image and a FDG uptake mean value is computed at each surface 
node of the medial surface (Fig. 7b, middle contour) in a 5 x 5 x 5 neighborhood 
(Fig. 7c). 

3 Displays of 3-D Functional Maps 

Fig. 8 a,b illustrate the FDG metabolic activity over the LV medial surface for 
2 cases. Right ventricular and epicardial surfaces are shown in transparency. 
Visualizations are performed using the VTK software library [9] . 
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(a) (b) 



Fig. 4. (a) 3-D prior biventricular model, (b) Corresponding mesh (1853 nodes, 7257 
elements) . 




Fig. 5. Model immersed in SA MR images. 




Fig. 6. Segmentation results of basal (BS) and mid- ventricular (MV) slice of (a) patient 
HI (BS), (b) patient HI (MV), (c) patient H2 (BS) and (d) patient H2 (MV). 
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(a) (b) (c) 



Fig. 7. (a) Intersection of endocardial (intrinsic contour), epicardial (outermost con- 
tour) and calculated medial (middle contour) surfaces with a MR SA image, (b) the 
same contours are shown in the corresponding registered PET SA image (medial con- 
tour in middle) and (c) the mean value of the FDG-uptake is computed at the medial 
surface nodes in a small neighborhood. 




(a) 



(b) 



Fig. 8. 3-D representation of the FDG PET uptake values of (a) patient HI and (b) 
patient H2. 



4 Discussion and Conclusion 

Using the described imaging protocol, we are able to combine MR anatomical 
and PET functional data for the same heart. The reliability of the result mainly 
depends on the accuracy of the PET-MRI registration. Visually satisfactory 
results were obtained with 9 over the 10 available cases. For one bad case, there 
was unexpected artifacts in FDG PET data. Two of the good biventricular 
deformable model results were illustrated in this paper. For some of the cases, 
slight manual intervention to move medial surface in PET image was needed. 
For better validation of the method phantom experiments or simulated images 
will be needed. 
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In [5], we have proposed a quantification of the registration accuracy based 
on a distance between segmented surfaces. Transaxial MR images were obtained 
with a snapshot technique during free respiration. In comparison, PET trans- 
mission images could be considered as an integral image over the acquisition 
time (about 10 minutes). This might lead to differences in the thorax shape 
and diaphragm position between these two imaging methods and causes error 
in this surface based registration method. Also twisting and contraction of the 
heart causes registration errors between gated cardiac MR and non-gated PET 
images. The registration of SA images was performed between end-diastolic cine 
MR images and ungated PET images. The end-diastolic time point of cine MR 
images was selected, because ungated PET images result from the uptake inte- 
gration over several heart cycle and mostly represent the diastolic phase. We are 
going to investigate the accurate validation of the registration approach through 
image simulations. 

In the segmentation procedure, the initialization step is crucial. We are cur- 
rently working on its automation. As the model is locally deformed, it has to 
be positioned quite close to the contours. This could be achieved by a rigid 
registration method (Iterative Closest Point [10] for example). 

One can also propose different alternative strategies to attribute PET-EDG 
uptake values to the model. This concerns our future work. We believe that 
the 3-D functional maps proposed in this paper are valuable tools to help the 
cardiologist in the interpretation of multi-modality cardiac imaging data for the 
better diagnosis of the myocardial viability. 
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Abstract. A 3D image processing method to increase reproducibility 
in the handling of PET myocardial perfusion studies is proposed. It is 
based on the 3D regularization of factor volumes, which are estimated 
by EAMIS (Eactor Analysis of Medical Image Sequences). The result- 
ing regularized factor volumes correspond to right and left cavities and 
to perfused tissues. They are then submitted to a C-means classifica- 
tion, looking for 4 clusters. Some rules of connectivity are applied to 
the resulting clusters to achieve the segmentation of heart cavities and 
myocardium. Thereafter, kinetics in the left ventricle and in the my- 
ocardium can be computed from the previously determined volumes of 
interest, and myocardial blood flow (MBE) is estimated by a conventional 
compart mental analysis. 

1 Introduction 

PET myocardial perfusion studies are primordial when absolute quantitation of 
myocardial blood flow is required [I]. Oxygen-15 labeled water (H2^^0) studies 
offer a large potential interest, since water is an ideal tracer to study perfusion. 
The short half-life of oxygen- 15 allows to repeat experiments in a short time, 
which is useful to study rest and stress conditions, for instance. However, the 
signal to noise ratio (SNR) of these studies is very low. That is the reason 
why some additional scans (using or ^^FDG) are needed to delineate the 

myocardium [2] . This has hindered the generalisation of such studies in a clinical 
context. To overcome this limitation, some a posteriori image processing methods 
have recently been proposed. They are mainly based on statistical data analysis 
[3], and more particularly on factor analysis methods, such as Factor Analysis 
of Medical Image Sequences (FAMIS)[4],[5] . In case of a first-pass myocardial 
perfusion study, FAMIS can synthesise a whole dynamic series into 3 or 4 kinetics 
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(called factors) and associated images (called factor images), which correspond 
to 1) right cavities, 2) possibly lungs, 3) left cavities and 4) myocardium and 
other perfused tissues. FAMIS algorithm consists first of estimating factors, then 
of computing factor images from raw data and factors, using a conventional least- 
squares method. The quality of factor images strongly depends on the noise level 
present in the raw data [6] . We have proposed a 2D regularisation of factor images 

[7] to improve the quality of factor images in case of raw studies with a low SNR. 
Its principle is to add some neighbourhood constraints in the reconstruction of 
the factor images, and so to improve the quality of factor images, especially 
regarding the SNR. The method has proved to be useful in the processing of 
H2^^0 myocardial studies [4]. The first applications were carried out on exams 
acquired with a multi-slice PET scanner. However the method does not take 
into account the features of up-to-date PET scanners, which provide data with 
a nearly 3D isotropic resolution. Moreover, with up to 63 slices, the data volume 
to be processed becomes huge and manual contouring of volumes of interest is a 
complex and fastidious task. 

This paper describes the developments that we propose to better process 
3D data. Eirstly, the generalization of the regularization method to 3D data is 
described. Secondly, a method for the automatic contouring of heart cavities and 
myocardium is proposed. The processing chain is illustrated on a clinical study. 
The potential interest of the method is discussed and possible applications in 
other modalities are proposed. 

2 Materials and Methods 

2.1 Materials: Dynamic Perfusion Studies 

Dynamic studies were acquired with an EC AT HR+ PET scanner, according 
to the protocol given in [2]. The injected dose was lowered to take into account 
the sensitivity of the PET scanner. The acquisition was performed in 2D mode, 
i.e. with the septa. Eor each time, volumes were reconstructed using a standard 
EBP algorithm. Pour dimensional data sets (three space dimensions plus time) 
were finally obtained. 

2.2 FAMIS Algorithm 

The EAMIS algorithm was extended to process 4D data sets, as described in 

[8] . Its consists in reducing the data set S{v,t), where v represents the voxel 

variable and t the time variable, to a limited number, of factors, Ck{t), and 
of associated factor volumes, according to the equation : 



Nk 

S{v,t) = '^Ck{t).Dk{v) + e{v,t). 

k=l 



( 1 ) 



The term e{v^ t) is the error term associated to each voxel at a specified time. 
The error is assumed to be zero-mean, and its covariance matrix is assumed to 
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be known. Positivity constraints are applied to both Ck{t) and terms, in 

order to be consistent with a physiological interpretation. 

Once the factors Ck{t) are estimated, the conventional estimation of the 
factor volumes Dk{v) is based on the minimisation of the term Ei : 

Nt 

t=l v=l 

Nt / Nk 

= EE ( -'^Ck{t).Dk{v) 

(=1 v=l V k=l 

Finally, the contribution Ck of the couple of factor and factor volume is 
computed according to the equation : 



( 2 ) 

(3) 






2.3 3D Regularisation of Factor Volumes 

The goal of the regularisation is to improve the quality of factor images, by 
introducing some a priori knowledge in their construction. The conventional 
way consists of setting some similarity constraints between neighbouring voxels, 
which has for mathematical formalism: 

N X) N XIV 

= E E haPfeK), Dk{v)), (5) 

v = l v*=l 

where N^v is the number of voxels, which are connected to the voxel v. Practi- 
cally, the six-connectivity is chosen: it is sufficient to introduce the neighbouring 
constraints and it enables us to reduce the number of computations which are 
required by the 26-connectivity. The he, function is the Hubert function: it is 
quadratic for small differences in image intensity and linear elsewhere. 

{ x^ for — a < X < a 

a{2x — a) for X > a 
—a{2x -h ol) for x < —a 

The regularised factor volumes Dj^{v) are estimated by minimising the term 
E: 

Nk ^ 

E = El ^ —E2k (6) 

k=l ^ 

where /i is a positive coefficient which sets the trade-off between the fit to the 
raw data and the a priori of smoothness in the factor volumes. 
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2.4 3D Classification and Segmentation of Heart Cavities and 
Myocardium 



A Nk component vector 



( 



Dijv*) D2 (v*) 

MaXv{Di(v)) ’ MaXv{D2(v)) ’ 



PN^jv*) 

Max-a(DNk(i^)) 



^ is associ- 



ated to each voxel i;*. These Ny vectors are submitted to a C-means algorithm. 
When Nk = 3, clusters are initialised by the four following vectors (1,0,0), 
(0, 1,0), (0, 0, 1), and (0.1, 0.1, 0.1). The resulting clusters are functional clusters 
and they do not match the anatomic regions. For instance, myocardial and hep- 
atic tissues are in the same cluster, since their functional behaviour is nearly the 
same in a perfusion study. To separate such functional structures, some further 
post-processing has to be applied to 3D clusters. The largest 3D connected vol- 
umes are first isolated. Then left cavities are separated from large vessels, such 
as the aorta. In a last step, a simple a priori geometrical knowledge is introduced 
to define myocardial tissue: myocardium is defined as tissue surrounding the left 
ventricle. A 3D spatial contiguity constraint is applied which ensures that the 
maximal distance between one myocardial voxel and one left ventricular voxel 
should be less than 4 voxels. 



3 Results 

For first-pass myocardial perfusion studies, FAMIS estimates generally three 
structures, which can be identified to right cavities, left cavities, and perfused 
tissues, which are mainly the myocardial and hepatic tissues in the field-of- 
view of the exam.. Figure 1 shows the regularised factor volumes which have 
been estimated using a = ^ _ q 25. The term E is minimised 

by the Flechter-Reeves’s algorithm. Convergence is generally obtained after 20 
iterations. 




Fig. 1. Regularised factor volumes (27 consecutive slices) estimated by FAMIS soft- 
ware. From left to right : right cavities, left cavities, and myocardium plus other per- 
fused tissues, such as liver. 
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The C-means classification identifies 4 clusters, which are defined as 1) right 
cavities, 2) left cavities, 3) myocardium and other perfused tissues, and 4) back- 
ground and lungs. The further processing of clusters which is based on geo- 
metrical considerations provides the segmentation which is shown by Figure 2. 
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Fig. 2. Results of the segmentation on the 27 consecutive slices. Surface view of seg- 
mented structures generated by the Anatomist software (SHFJ, CEA). Right cavities 
are in blue colour, left cavities in red colour and myocardium is in light maroon colour. 



4 Discussion 

The interest of the regularisation of factor images has already been demonstrated 
[7]. The generalisation of the method to 3D data is necessary to take into ac- 
count the specific 3D isotropic acquisitions realized by new PET scanners. The 
classification is a way to segment the different cardiac structures based on the 
regularised factor images or factor volumes. The C-means algorithm is robust 
when an adequate initialisation is performed. In our case, this initialisation can 
be easily performed, since the clusters can be a priori defined : right cavities and 
left cavities, myocardium, and background. In some cases, it may be useful to 
add a fifth cluster, corresponding to lungs and which can be initialised by the 
(0, 0.3, 0.3) vector. Compared to threshold-based methods, which need a priori 
calibration, our method is more simple to adapt to different cameras. Another 
interest of the method is that it takes into account the functional information 
which is in the dynamic first-pass data. As a consequence, right cavities, left 
cavities and perfused tissues can be isolated during the same procedure. The 
knowledge of the left ventricular position is useful to segment the myocardium. 
Thereafter, the geometrical constraint which has been introduced to position 
the myocardium is sufficient to separate the myocardium from the liver. The 
potential applications of this segmentation are not restricted to H 2^^0 studies. 
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It could be applied to PET studies, which are largely used for assessing 

myocardial perfusion. Moreover, some applications of this principle have been 
tested on MR first-pass perfusion studies. The goal of the proposed segmentation 
is not to solve the partial volume effect : some more precise a priori information 
would then be necessary to take it into account. Our recommandation is that the 
contamination of myocardium by the left cavity has to be taken into account by 
the compartmental model, as it is done when manual segmentation is performed. 
The limits of the algorithm in case of the different pathologies, which introduce 
a reduction of MBF, should be defined. 

5 Conclusion 

An automatic method to segment first-pass myocardial perfusion studies has 
been proposed. It takes advantage of the presence in the same study of heart 
cavities and myocardium. Although the different structures are not well sepa- 
rated on raw data, the processing based on factor analysis and followed by a 
classification allows us to separate them. The method needs no a priori calibra- 
tion and its principle can be applied to a large variety of first-pass myocardial 
studies. 
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Abstract. Aim of the study was to test the feasibility of cine-NMR for assessment of 
the infarcted rat heart and to compare the results to established methods. Thereafter, 
the value of cine-NMR was tested in studies investigating interventions to change 
the course of the remodeling process. NMR was performed for determination of left 
ventricular (LV) volumes and mass, Ml-size and cardiac output. After MRI rats 
underwent conventional hemodynamic measurements for determination of cardiac 
output and LV volumes by electromagnetic flowmeter and pressure-volume curves. 
LV wet weight was determined. Ml-size was determined by histology. MRI- 
acquired Ml-size (18.5±2%) was smaller than histology (22.8±2.5%, p<0.05) with 
close correlation (r=0.97). There was agreement in LV mass between MRI and wet 
weight (r=0.97, p<0.05) and in the MRI- and flowmeter measurements of cardiac 
output (r=0.80, p<0.05). Volume by MRI differed from pressure-volume curves with 
good correlation (r=0.96, p<0.05). In conclusion, cine-NMR is a valuable diagnostic 
tool applicable to the rat model of MI. Being non-invasive and exact it offers new 
insights in the remodeling process after MI because serial measurements are 
possible. 



1 Introduction 

The rat heart infarct model has proved to be predictive for cardiac remodeling after 
myocardial infarction (MI) in patients and for therapy [1]. Being noninvasive and 
exact nuclear magnetic resonance (NMR) can offer new insights allowing serial 
investigations of rat cardiac morphology and function. Aim of the study was a 
comparison of data acquired by NMR and conventional, established methods [2, 3]. 
After completion of the validation experiments the technique was applied to several 
interventional studies in the rat heart infarction model, as impact of transmyocardial 
laser revascularisation on remodeling after MI, influence of androgene levels on post 
infarct remodeling and influence of a cerivastatin therapy on left ventricular 
remodeling. 



2 Methods 

In group A 8 controls and 8 rats 16 weeks after myocardial infarction (MI) were 
investigated on a 7 T-Biospec (Bruker, Germany) using an ECG-triggered Cine- 
LLASH-sequence. A home built rat size whole body coil was used as transmitter and 
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a surface coil as receiver. For imaging, we used an ECG-triggered fast gradient echo 
(FLASH) cine sequence [ 4]. Flip angle was 30 to 40°, echo time was 1.1 ms, 
repetition time 3.2 ms. 12 frames per heart cycle were obtained. The total acquisition 
time (TAT) for one cine sequence was 40 to 50 s depending on heart rate (TAT=128 
phase encoding steps * 2 averaging steps * length of one heart cycle). The 
measurements were averaged four times to increase signal to noise ratio. 15 to 18 
contiguous ventricular short axis slices of 1 mm thickness were acquired to cover the 
entire heart. Total scan time was in the range of 15 min. With a field of view of 3 to 4 
cm and an image matrix of 128 by 128 in plane resolution was 230 to 310 pm. Data 
analysis was done using an operator- interactive threshold technique. Myocardial and 
ventricular slice volumes were determined from end-diastolic and end-systolic images 
by multiplication of compartment area and slice thickness (1mm). Total volumes were 
calculated as sum of all slices. Left ventricular mass was calculated as LV myocardial 
volume multiplied by the myocardial specific gravity (1.05 g/cm^). Myocardial infarct 
size (Ml-size) was determined for every slice as the myocardial portion with 
significant thinning and akinesia or dyskinesia during systole. Ml-size was calculated 
by dividing the sum of the endocardial and epicardial circumferences occupied by the 
infarct by the sum of the total epicardial and endocardial circumferences of the left 
ventricle. This method was adapted from the histological determination of Ml-size 
first described by Pfeffer et al. [2]. 

After NMR hemodynamic measurements were performed for determination of the 
same parameters by electromagnetic flowmeter and postmortal pressure -volume 
curves (simultaneous infusion at 0.76ml/min and pressure recording, end-diastolic 
volume is infused volume at in vivo measured end-diastolic pressure) [2, 3]. The left 
ventricles were weighted. In group B NMR-acquired and histological Ml-size (picric 
red dye after formalin fixation) [2] were compared in 26 rats. 



3 Results 



As shown in the table, good correlation between NMR and invasive techniques was 
found. Significant differences between methods were detected in infarct size, end- 
diastolic volume and ejection fraction. 





NMR 




conventional method 


r 


Ml-size % 


18.5±2 




22.8±2.5='^ 




1-0.97 




MI 


control 


MI 


control 




LV mass (mg) 


865.1±39.2 


537.6±19.6 


865.1±41.3 


540.3±18.4 


1-0.97 


EDV (III) 


737.0±70.5 


343.9±8.4 


671.1±64.1* 


262.7±12.8* 


r=0.96 


SV ( III) 


259.1±10.4 


224.4±7.8 


245.8±13.5 


221.2±9.2 


r=0.70 


EF (%) 


36.7±2.4 


65.3±1.6 


41.9±3.8* 


84.4±2.3* 


r=0.97 



Data is given in meantSLM. r: correlation coefficient of methods (MI plus control), 

^ p < 0.05 vs. NMR of the same group, LDV: end-diastolic volume, SV: stroke 
volume. 
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5 Conclusion 

The study showed good agreement of NMR-acquired data with established methods 
for LV mass and stroke volume. With a difference of 5% both methods of Ml-size 
determination correlate well. The difference to the NMR-method is due to 
inhomogeneity of shrinkage caused by formalin fixation of the excised hearts. There 
was good correlation of postmortal determined to in vivo measured end-diastolic 
volume. The significant difference was bigger in the controls than in Ml-rats ( 23.4% 
and 9.0%) and is attributed to the absence of enregy dependent, active relaxation in 
hearts arrested by potassium chloride. In conclusion, application of NMR to the rat 
model of myocardial infarction is promising to shed new light on the 
pathomechanisms of cardiac remodeling, because changes during this process can be 
monitored in intact individual animals. Serial characterisation of infarct size allows 
investigation of ongoing changes of the Ml-scar. 




Fig.l. Diastolic long axis frame of control rat. This is one of twelve frames covering the cardiac 
cycle (average length in the rat 300 ms), which can be animated as a movie. Image shows left 
ventricle with inflow and outflow tract. 




Fig.2. Systolic short axis frame of rat 12 weeks after anterolateral myocardial infarction. 
Infarcted region has lost its function and is characterized by a thinned wall. 
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4 Application Studies 

4.1 Transmyocardial Laser Revascularisation Improves Perfusion but 

Enhances Left Ventricular Remodeling in Rats after Myocardial Infarction 

Hypertrophy following MI causes diffuse hypoperfusion of the remote myocardium. 
We therefore used the MI rat model for assessment of transmyocardial laser 
revascularisation (TMLR). 8 weeks after coronary artery ligation rats had Cine-NMR, 
thereafter TMLR was done. A selfdesigned HolmiumiYag laser system emitting a 
wavelength at 2.1 pm (pulse energy 4J, pulse duration 1 ms, 365|um fiber diameter) 
was placed in soft contact to the remote myocardium. 4 weeks after TMLR a second 
Cine-NMR including dobutamine-stress (10 pl/kg/min via tail vein) was follwed by 
high resolution spin labeling perfusion imaging (pixel size 140 pm). For assessment 
of perfusion changes induced by TMLR treatment, spin-labeling MRI was performed 
as described previously [5] in a 11.75 Tesla wide bore magnet (AMX 500, Bruker, 
Karlsruhe, FRG). The principle of the measurement is that spins of a selected imaging 
slice are inverted and Tj relaxation of these spins towards equilibrium is observed. 
Due to perfusion non-inverted (equilibrium) spins flow into the slice. This results in a 
decrease of the apparent proton relaxation time in tissue. Perfusion imaging was 
performed in isolated retrogradely perfused hearts obtained from rats of the TMLR 
group (n=10). 

TMLR-areas were located by corresponding H&E histology. TLMR-areas were 
better perfused (3.89±0.83ml/min/g at rest and 2.29±1.06ml/min/g at nitro stress, p < 
0.05 both). Wall thickening at rest was higher in TMLR regions than in control 
regions of the same heart. (57+7 vs 37+6%, p<0.05). 





Control 


TMR 


Delta EDV (|al) 


24.6±16.7 


81.7+15.7^ 


Delta LV mass (mg) 


54,5+19.2 


124.1+30.7^ 


EF at 12 weeks (%) 


40+2 


38+2 


Cl at 1 2 weeks (ml/kg/min) 


232.2+12.8 


242.4+12.8 


EF at 1 2 weeks stress (%) 


43+5 


54+5^ 



Mean+SEM. =^p < 0.05 vs control, delta: difference 8 vs 12 weeks after Ml. EDV: 
end-diastolic volume, LV : left ventricle, EF: ejection fraction, Cl: cardiac index. MI- 
size 27+2% both groups. 

Using high resolution perfusion NMR, improvement of perfusion by TMLR was 
visualized. Due to the retrograde perfusion setting of the isolated heart this 
improvement cannot be caused by flow from inside the LV cavity via open channels. 
Translation of the results into the patients situation is not straight forward, however, 
some conclusions can be made. One explanation for the discrepancy of improved 
angina in patients after TMLR, but the failure of major studies to show improved 
perfusion may be insufficient spatial resolution of the applied imaging techniques. 
Addititonally, patients with impaired microcirculation may be successfully treated 
with TMLR. In relation to the left ventricular size in the rat the diameter of the laser 
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fiber was rather large, however, the finding of enhanced LV remodeling should lead 
to additonal caution for the application of TMLR in the situation of post MI 
remodeling. 



4.2 Influence of Testosterone on Cardiac Remodeling after Myocardial 
Infarction in Rats 

Several facts suggest, that left ventricular remodeling after myocardial infarction may 
be influenced by androgen levels [6-8]. There is evidence of testosterone receptors in 
the heart. Further, it has been shown that testosterone use in bodybilders leads to 
increased left ventricular mass and rhythm abnormalities. Additionally, there is a 
large body of evidence for gender differences in pathophysiological processes in 
cardiovascular disease, for instance a lower incidence of myocardial infarction in 
premenopausal women. Aim of the study was to assess the influence of androgen 
levels on left ventricular remodeling after myocardial infarction (MI) in the intact rat. 
Wistar rats were assigned to 3 groups: surgical removal of testes 2 weeks prior to MI 
(ORX), Placebo (PL) and Testosterone-undekanoat-Injection (TUD, 500mg 
intramuscular 2 weeks prior to and 2 weeks after MI). NMR-measurements were done 
2 and 8 weeks after left coronary artery ligation. Left ventricular (LV) mass and 
volumes, cardiac output (CO), ejection fraction (EF), wall thickness (WT) and - 
thickening (SWT) were determined. In addition, a transversal gradient echo of the 
paravertebral region was acquired to evaluate the growth of sceletal muscle. 





ORX (n=13) 


Placebo (n=l 3) 


TUD (n=9) 


MI size 2 weeks (%) 


36.4±2.0 


38.8±2.9 


37.7±1.8 


A EDV (pi) 


+276.0+33.2 


+263.9+41.2 


+245.6+51.2 


A LV mass (mg) 


+82.7±29 


+170.7±43.6 


+274.2±65.9* 


CO 8 weeks (ml/min) 


104,8±7.2 


103.5±10.3 


102.0+4.7 


EF 8 weeks (%) 


36+3 


33+4 


35+2 


SWT 8 weeks (%) 


66±7 


34±3^^ 


30±4* 


PM (mm^) 


145.3±4.4 


168.0±4.5^ 


187.8+5.3*t 



Mean±SEM. p < 0.05 vs. ORX, t 4 p < 0.05 vs. Placebo. A: difference from 2 to 8 
weeks after ML EDV: end-diastolic volume, PM: area of parevertebral sceletal 
muscle. 

Dilation and global LV function after MI was not changed by manipulation of 
androgen levels, but hypertrophy was. Regional wall thickening of the remote region 
was preserved in ORX. In ORX growth of sceletal muscle was impaired and 
accelerated in TUD. 



4.3 Impact of Cerivastatin on Cardiac Remodeling after Myocardial Infarction 

Statin therapy has been shown to lower the incidence of myocardial infarction via 
positive effects in atherosclerosis [ 9]. Recently, it has been reported that statins have 
effects other than lipid lowering, and treatment after MI may be beneficial not only 
due to secondary prevention of reinfarction [10]. The influence of long-term 
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cerivastatin therapy on left ventricular remodeling was assessed in rats with chronic 
myocardial infarction. In addition, the contribution of endogenous nitric oxide (NO) 
formation was assessed in animals treated with the NO synthase inhibitor L-NAME in 
combination with hydralazine (to avoid the increase in blood pressure induced by L- 
NAME). Using NMR it was possible to serially determine left ventricular mass as the 
most exact parameter of hypertrophy . 

NMR-measurements were done 4 and 12 weeks after left coronary artery ligation. 
Rats were treated either with placebo, with cerivastatin (C, 0.6 mg/kg body weight) 
starting on the 7^^ postoperative day as dietary supplement, or cerivastatin plus L- 
NAME (76 mg/ 100 ml) and hydralazine (8 mg/ 100 ml) in drinking water (CLH)) 
daily for 12 weeks. 





Sham/ 

Placebo 

(n=8) 


MI/ 

Placebo 

(n=8) 


Cerivastatin 

(n=ll) 


CLH 

(11=11) 


MI size 4 weeks (%) 


0 


31.512.7 


31.911.4 


31.912.8 


A EDV (pi) 


+26.8±20.4 


+108.7128.8^ 


+126.6120.5^ 


+ 173.7125.U 


A LV mass (mg) 


+25.318.6 


+235.3133.7^ 


+59.8120.5t 


+239.5116.0^ 


CO 12 weeks 


73.013.6 


76.112.9 


95.814.8^1* 


69.312.8:1: 


(ml/min) 

EE 12 weeks (%) 


70.112.0 


36.213. U 


40.713. U 


33.012.6^$ 


A EDWT (mm) 


0.010.06 


+0.510.06^^ 


-0.0810.06t 


+0.410.06*1: 


SWT 12 weeks (%) 


80.916.9 


32.514.4^ 


48.414.5t 


2912*1: 



Mean±SEM. p < 0.01 vs. Sham/Placebo, t P < 0.01 vs. Ml/Placebo, f p < 0.01 vs. 
Cerivastatin. A: difference from 4 to 12 weeks after ML EDV: end-diastolic volume. 



Administration of Cerivastatin attenuated left ventricular hypertrophy after MI. 
Dilation was not changed. Cardiac output and wall thickening were higher in 
Cerivastatin treated rats. The positive effects of Cerivastatin were completely 
abolished by NO synthase-inhibition. 



4 Discussion 

In the validation study NMR was shown to correlate closely to invasive or ex vivo 
techniques that are typically used in the rat model of myocardial infarction. With the 
advantage of its non-invasive character, NMR allows for serial investigations of the 
same animal. Hence, changes due to therapy or interventions can be monitored 
closely and more exact. Small differences can be detected, and the number of animals 
to be sacrified for research purposes can be reduced. 

In conclusion we could show that Cine-MRI is a valuable diagnostic tool 
applicable to the rat model of myocardial infarction, which is in good agreement with 
existing analytical methods. It allows for reliable and reproducible assessment of 
cardiac morphology and function making it the method of choice for volumetric 
quantification. In contrast to echocardiography NMR-volumetry works without 
geometrical assumptions and therefore offers a high degree of accuracy in 
measurements of cardiac output and LV volumes; this is especially of importance in 
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hearts deformed by asymmetric dilation after ML NMR can be applied to assess the 
effects of therapeutical or other agents on infarct size, cardiac geometry and function 
with high precision. Thus NMR bears a great potential to substantially contribute to 
the understanding of the underlying pathomechanisms in the development of chronic 
heart failure and can help to evaluate new therapy options. 
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Abstract. A modular data fusion system based on the Dempster-Shafer 
framework is presented. This system allows the building of any architec- 
ture of fusion by chaining elementary modules. Two types of modules are 
used, the numerical to symbolic conversion modules and the combination 
modules using logical rules of combination. The uncertainty is modeled 
by a basic belief assignment or the plausibility of each hypothesis. We 
applied our system to assess the Left Ventricular (LV) myocardial viabil- 
ity. The parameters taken into account are the LV contractile function 
extracted from tagged Magnetic Resonance Images (MRI) and the glu- 
cose metabolism rate obtained by Positron Emission Tomography (PET) 
imaging. The variables of interest are defined by the medical experts. The 
results are displayed on polar maps to give a geometrical information of 
the potential lesions. 



1 The Modular Fusion System 

Medical diagnosis such as myocardial viability results from the integration of 
many complex phenomena which can not be precisely described by the special- 
ists. Moreover, medical knowledge continuously evolves and research results lead 
to take into account new parameters. Our modular fusion system is designed to 
be adaptive to this uncertainty and this evolution. It allows medical specialists to 
experiment and analyze new architecture of fusion by combining the parameters 
according to their choice. This system enables to build any architecture of fu- 
sion. An architecture of fusion defines the scheme of data processing. It consists 
of a set of chained modules, classified in two types : the numerical to symbolic 
conversion modules (named Num2Sym modules) and the combination modules. 
The first modules are used for symbolic interpretation of the numerical reports 
from the sensors. The second ones perform the combination of these symbolic 
data in order to get synthetical information useful for diagnosis. 

Running the algorithms of these modules needs a priori knowledge which 
is : (a) in case of Num2Sym modules, the model of the symbols defined on the 
reference of the numerical reports, (b) in case of combination modules, the logical 
rules of combination which define the symbolic values of the output from the 
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symbolic values of the inputs. In this work, a priori knowledge will be provided 
by medical experts. In figure 1, an example of a simple architecture of fusion 
is presented for two inputs parameters X and Y . The numerical inputs are 
represented in italic in grey ellipsoids. These inputs are converted to symbolic 
data by Num2Sym modules, giving the symbolic data in white ellipsoids. These 
data are finally fused to get the output parameter Z. 




Fig. 1. Simple architecture of fusion 



2 Dempster Shafer Framework 

Our system is based on the Dempster-Shafer framework [1,2,3], well-adapted to 
model the uncertainty of the data and the rules of decision. 



Model of the numerical to symbolic conversion. Each numerical pa- 
rameter X taking x value in the interval Ix is linked to a set of hypothesis 
Qx = {Hx/i = 1, N the number of hypothesis. It is proposed to define a 
trapezoidal function assigned to each hypothesis on the interval Jx, 
but also to each proposition stemmed from the logical union of several hy- 
pothesis. In fact, the medical expert always defines ordered hypothesis on the 
interval Ix with inaccurate thresholds, so we propose to model the doubt be- 
tween to consecutive hypothesis by a composed hypothesis or proposition such 
as where U corresponds to the union of the two hypothesis, such as 

represented Figure 2. The function called basic belief assignment, models the 
confidence that a numerical value x £ Ix belongs to the different propositions 
defined by subsets of i?x- The definition of rrix is : 

m, : 2^- ^ [0, 1] 

A ^ = /a^{x) 

We deliberately choose a simple model, because it is well-adapted to translate 
the medical expertise. In our application, propositions given by experts are or- 
dered on the numerical space of definition Ix, only the thresholds between two 
consecutive propositions are not accurate. 



Model of combination rules. The combination of two symbolic variables 
gives a third symbolic variable that must be interpreted by the medical experts. 
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Fig. 2. Model of the numerical - symbolic transformation 



The basic combination rules are extensively described by a decision table from 
the a priori knowledge of the clinical experts. Let f2x = f2y = 

{Hy^ Hy: Hy} aud Qz = {Hzi be the sets of hypothesis respectively 

related to the symbolic variables X, Y and Z. Table 1 exhibits an example of 
rules determining the output variable Z from the combination of X and Y. 



Table 1. Basic rules of combination defined from X and Y for the determination of Z 



X hypothesis 





Z 




Hi 


Hi 


Y 


Hi. 


Hi 


Hi 


Hi 


hypothesis 


Hi 


Hi 


Hi 


Hi 




Hi 


Hi 


Hi 


Hi 



The combination is based on the unnormalised form of the Dempster’s con- 
junctive rule. This rule allows the combination of two distributions defined on 
the same space of discernment. Nevertheless, it is possible to extend this rule, if 
a logical table of combination such as Table 1 is known. An illustration is given 
with the following simple example. Let 

m, : 2 ^- ^ [ 0 , 1 ] 

my : ^ [0, 1] 

be the basic believe assignments of the numerical parameters X and Y. rrix and 
rriy can be redefined by and m'y on the discernment space 

nx®nY = {{Hj, n h^), {Hj, n h^), {h\ n h^), {h\ n hD) 

n H^) U {H\ n H^)) 

m,{H\ U H\) = n U {H\ n H^) U {H\ n H\.) U {H\ n H^)) 
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The same extension can be made with rriy to m'y. Thus, m'^ and m'y are 
defined on the same space of discernment. The Dempster’s combination rule can 
be used given m'^ = 0 where m'^ is defined on for the two 

measures x and y. 

We can notice that this extension increases the dimension of the space of 
discernment. Dim{Qx) = 2 and Dim{QY) = 2, so Dim{Qx^^Y) = 4. In many 
cases, the expert gives an interpretation of the resulting hypothesis, defining the 
hypothesis of the result variable Z. For instance. 

This interpretation can be described by a logical array as shown in Table2. 
The transformation of defined on Qx 0 to mz defined on Qz is made 
by using equations (4). 



Table 2. Interpretation of the resulting hypothesis 



X hypothesis 





Z 


Hk 


Hi 


Y 




Hh 


Hi 


hypothesis 


Hi 


Hi 


Hi 



Practically, it is sufficient to extend the Table 2 from i? to 2^ to get the new 
Table 3 that represents the Tab function, and to transform the Dempster’s rule 
to: 

mz{A) = y] m^{B).my{C) 

A=Tab{B,C) 



Table 3. Extension of the resulting hypothesis 



X hypothesis 





Z 


Hi 


Hi 


HluHl 


Y 


Hi 


Hi 


Hi 


HluHl 


hypothesis 


Hi 


Hi 


Hi 


HluHl 




HluHl 


HluHl 


HluHl 


Hi U Hi U Hi 



Reports processing. Both models developed in the previous section are com- 
pletely independent of the reports of a particular patient. They just model the 
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knowledge of the medical experts. The reports processing uses the previous mod- 
els. Given a numerical value x of the variable X, a basic belief assignment rrix 
is computed from the model of conversion previously defined in the module 
Num2Sym associated to this variable. The combination step needs the defini- 
tion of a combination table. The belief assignment ruz of the output variable Z 
is computed from the Dempster’s conjunctive rule in unnormalized form. 

3 Application to the Diagnosis of Myocardial Viability 

Myocardial viability. Myocardial viability assessment is an important field 
of research in cardiac diseases. Viable tissue is under-perfused tissue that can 
recover after blood flow reestablishment. The quantification of the viable tissue is 
of major importance in medical diagnosis since it helps to decide whether or not 
a patient will benefit from a revascularisation procedure [4] . In order to make an 
automatic accurate viability analysis, complementary diagnosis parameters have 
to be integrated in to the system. It has been shown that the combination of the 
glucose metabolism and the LV contractile function gives a good assessment of 
the viability [5,6]. Contractile function analysis is obtained from tagged Magnetic 
Resonance Imaging (MRI). Glucose metabolism rate is provided by a Positron 
Emission Tomography (PET) scanner. 



Input parameters. MR images of the heart are acquired using a Siemens Vi- 
sion 1.5T scanner with 35ms temporal resolution (8mm slice thickness, 45 and 
135 degree tag orientations, EOV 280mm). The investigation of LV function in 
tagged MR images is performed with a cardiac analysis package {FindTags, 
John Hokkins University, Baltimore). Two parameters are considered for the 
assessment of the LV contractile function: the Maximum circumferential Short- 
ening (MS) of deformation tensor and the Angle (A) formed by displacement 
vector versus the circumferential direction. Both parameters are computed on 
rest images {MSr, Ar) and stress images {MSs, As) at end diastole. PET vol- 
ume (64 slices) are acquired on a HR+ tomograph (Siemens Erlangen) in 3D 
mode with intrinsic resolution of (4mm x 4mm x 4mm). One hour after initia- 
tion of euglycemic hyperinsulinemic clamp infusion, 3 to 4 mGi of 18EDG are 
injected as a slow bolus. The studied volume is taken 45 min after this infusion 
with 15mn acquisition time. The uptake rate (FDG) (glucose metabolism pa- 
rameter) of any myocardial region is obtained by counting procedure. A major 
problem in fusion system is data registration. In our application, this step is 
previously performed using a semi-automatic software {Cardiofuz^ GREATIS, 
Lyon). Our modular fusion system will use five input parameters: MSr, Ar, 
MSs, As and FDG. The format of these data will be uniformed and fits well to 
a standard polar map display. 



Architecture of the fusion system. Such as presented Eigure 3, the numer- 
ical inputs of MSr, Ar, MSs, As and FDG are converted into symbolic data 




Dempster Shafer Approach for High Level Data Fusion 109 



by Num2Sym modules, giving the symbolic data in white ellipsoids. These data 
are finally fused two by two until to get the needed results which is the viability 
(Viah). Intermediate variables are derived from the initial ones: the contractile 
function at rest {CRr), the contractile function under stress {CFs), the inotropic 
function based on CFr and CFs comparison (IF) and the viability Viab. 




Fig. 3. Architecture of fusion 



Models of the numerical to symbolic conversion. Num2Sym modules 
are configured in order to make the conversion of the numerical values of MSs, 
MSr, Ar, As and FDG into symbolic data. The hypothesis and the basic belief 
assignments of each numerical variables are defined with medical experts. The 
symbolic variables obtained from the numerical variables are described in Table 
4. For instance, the MSs variable is associated by the medical expert to three 
hypothesis, then Qmss = Each trapezoidal function /^v, /z, fp and 

of course /atuz, fzup is also defined by the expert. These definitions allow the 
system to make the numerical/symbolic conversion as presented section 2. 



Table 4. Single hypothesis of the symbolic input variables 



MSs 


As 


MSr 


Ar 


FDG 


N 


Negative 


N 


Negative 


S 


Small 


S 


Small 


S 


Small 


Z 


Zero 


Z 


Zero 


L 


Large 


L 


Large 


M 


Medium 


P 


Positive 


P 


Positive 










L 


Large 



Models of combination. Table 5 describes the symbolic variables CFr, CFs, 
IF and Viab derived from the symbolic input variables MSr, MSs, Ar, As and 
FDG. The basic rules of combination are represented by logical arrays. Table 
6 summarizes these combination rules for all symbolic output variables used in 
the architecture of fusion. 
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Table 5. Single hypothesis of the symbolic output variables 



CFs 


CFr 


IF 


Lr 


Lower 


ADK 


Akinetic 


ANI 


Akinetic at rest and no improvement 


Eq 


Equal 


HK 


Hypokinetic 


HNI 


Hypokinetic at rest and no improvement 


Hr 


Higher 


NL 


Normal 


FI 


Function improvement 










I 


Ischemic 










R 


Remote 



Viab 


Er 


Error : EDG /IE Mismatch 


V 


Viable 


NE 


Necrosis 


I 


Ischemic 


Md 


Maimed 


R 


Remote 


MV 


Metabolic Viable 







Table 6. Basic rules of combination used by the architecture of fusion 



MSr 



MSs 



CFr 



IF 


ADK 


HK 


NL 


Lr 


ANI 


HNI 


I 


Eq 


ANI 


HNI 


R 


Hr 


El 


El 


R 



CFs 


N 


Z 


P 


N 


Hr 


Hr 


Hr 


Z 


Lr 


Eq 


Hr 


P 


Lr 


Lr 


Lr 



CFr 


S 


L 


S 


HK 


NL 


L 


ADK 


ADK 



FDG 



IF 



Viab 


ANI 


HNI 


FI 


I 


R 


S 


Ne 


Md 


Er 


Er 


Er 


M 


Md 


Md 


V 


I 


R 


L 


MV 


MV 


V 


I 


R 



For instance, the belief basic assignment of the MSs and Ar variables are 
combined using the method presented in section 2 to give belief basic assignment 
of CFs variable. The space of discernment of CFs is then f2cFs = {Lr, Eq^ Hr}. 

Results. The medical experts need to know where a risk may occur, so we chose 
to show them the plausibility of each hypothesis of the final variable which is the 
viability. The plausibilities of each hypothesis presented in Table 5 are computed. 
The results obtained from the data of a particular patient are displayed on polar 
maps Figure 4. 

4 Conclusion 



A modular data fusion system was developed with Dempster- Shafer framework. 
It is based on modules of two types: the Num2Sym modules parameterized by 
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Fig. 4. The polar maps of the plausibility for the different hypothesis of the viability 




the thresholds and the frame of the conversion and the Combination modules 
easily definable to fit to the specialists knowledge. Thus, any architecture of 
fusion may be built from these modules. 

This system is a sort of rules base system, but taking into account the uncer- 
tainty on the data, modeled by a basic belief assignment. Similar results might 
be obtained with a fuzzy rule base system. The choice of the Dempster Shafer 
theory has been approved by the medical expert because of its ability to model 
the doubt. 

An architecture of fusion has been implemented to assess the myocardial 
viability. The diagnosis results from the processing of five input parameters CFr, 
CFs, Ar, As and FDG associated to LV contractile function and the glucose 
metabolism. It is displayed on polar maps, using the mass of evidence or the 
plausibility of each hypothesis. In a future work, we will assess the accuracy of 
the diagnosis by comparing our results to the real measurements performed on 
30 patients before and after a revascularisation procedure. 
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1 Introduction 



Perturbations in ventricular mechanical loading can be arrhythmogenic and have been 
associated with sudden cardiac death in patients suffering from congestive heart 
failure, dilated cardiomyopathy, or ventricular volume overload (1-3). Stretch- induced 
changes in action potential propagation or repolarization could provide a mechanism 
for mechanically induced arrhythmias. However, there is a paucity of information 
regarding the effects of altered load on conduction velocity. The few existing reports 
present an unclear picture; some of the discrepancies may be due to the varying 
techniques used. 

Stretch activated channels (SACs) have been identified as a potential cellular 
mechanism for mechanoelectric feedback (4). Mechano-sensitive ion channels have 
been identified in the ventricular myocardium of several species, both non-specific 
cation or specific channels (5-7). The aminoglycoside antibiotic streptomycin has 
been reported to block stretch activated channels (8). 

Past studies of cardiac mechanoelectric feedback have all used contact electrodes 
to measure the electrical activity of the heart. The required physical contact between 
the measuring device and the myocardium leaves the measurements open to the 
possibility of mechanically induced electrical artifact (9). Optical mapping provides a 
non-contact means of measuring cardiac electrical activity through the use of a 
voltage-sensitive fluorescent dye. The spatial resolution from this technique is higher 
than that of conventional electrode arrays and is thus particularly desirable for 
conduction velocity measurements. In the present study, we employed optical 
mapping to investigate the effects of increased left ventricular loading on apparent 
epicardial conduction velocity and action potential duration. Streptomycin was used 
to assess if any of the observed electrophysiological changes might be governed by 
stretch activated channels. A preliminary report of this study has appeared in the 
proceedings of the 2001 ASME Summer Bioengineering meeting (10). 

Computational models are also valuable tools for investigating cardiac 
electromechanical interactions. By adding stretch activated currents to ionic models of 
the cardiac myocyte action potential, the integrated effects on stretch on cardiac 
electrophysiology can be investigated (11). These models can then be integrated into 
three-dimensional continuum analyses of ventricular mechanics and action potential 
propagation. This work has all been previously reported in full or preliminary reports 
(10, 12-17). 
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2 Experimental Methods 

Experiments were conducted in isolated Langendorff-perfused rabbit hearts. A 
balloon in the left ventricle was connected to a fluid-filled pressure transducer and a 
volume infusion pump. Following aortic cannulation, hearts were perfused with 
Tyrode’s solution and allowed to actively contract. The initial volume in the balloon 
was adjusted to achieve an end diastolic pressure (EDP) of ~0 mmHg. Hearts were 
perfused with a bolus of the voltage-sensitive dye, DI-4-ANEPPS and loaded from 0- 
30 mm Hg left ventricular pressure at a constant rate. Just prior to the acquisition of 
optical measurements, the electromechanical uncoupler 2,3 butanedione monoxime 
(BDM, 12.5 mM) was also added to the perfusate to prevent motion artifact. The heart 
was paced at a cycle length of 300 msec from the left ventricular epicardium near the 
apex. The loading protocol was repeated after the perfusate was switched to one 
containing 200 ‘M streptomycin. 

The hardware configuration for optical mapping has been described previously 
(18). Fluorescence images of the LV free wall (lateral view) were captured with a 
digital CCD camera (Dalsa, Waterloo, Ontario) at a speed of 399 frames per second 
and a resolution of 128x128 pixels. The image processing and analysis of the optical 
data has been detailed elsewhere (18, 19). After the signals were filtered, activation 
times were identified as the time at the maximum first derivative (dF/dt)max of the 
optical action potential upstroke. Time of repolarization was computed by 
determining the time at which the optical action potential had recovered 20% and 
80% from its peak. The difference between the repolarization and activation times 
was taken to be the action potential duration (APD). 

The epicardial geometry of each heart was reconstructed from orthogonal biplane 
images of the lateral and posterior views, and a finite-element surface was then fitted 
to the boundary data yielding a three-dimensional model of the epicardium. The 
activation times, maximum derivatives of the upstroke, and action potential duration 
times were then mapped on to the surface of the finite element model and fit by linear 
least squares. This allowed the electrophysiological variables to be expressed as 
functions of the three-dimensional coordinates of the epicardial geometry. 



3 Results 

Epicardial activation times were reversibly prolonged during loading (Fig. 1). Total 
mean apparent conduction velocity averaged over all distances decreased from about 
400 mm/sec in the unloaded pre-control state to less than 300 mm/sec in the loaded 
state (p < 0.05). By ANOVA, the interaction effect between load state and distance 
from pacing site was not significant suggesting that the decrease in conduction 
velocity was not changed at the different locations on the surface of the heart. There 
was no significant change in apparent conduction velocity between the two unloaded 
states. (dF/dt)max of the optical action potential upstroke decreased about 5%, but a 
decrease of this amount or greater would be expected simply as a result of slowed 
conduction, so decreased membrane excitability is unlikely to explain the observed 
stretch- dependent slowing. Adding streptomycin to the perfusate did not change the 
overall conduction velocities or the response of apparent conduction velocity to load. 
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suggesting that SACs may not play a primary role in this effect. There was no 
significant difference in the relative decrease in apparent conduction velocity with 
respect to distance from the pacing site, suggesting that the effect may be present 
transmurally throughout the myocardium. 




Before 

Loading 




During After 

Loading Loading 



Fig. 1. Epicardial activation times in unloaded (left and right) and loaded (center) rahhit left 
ventricle. From (12) 

Total mean APD 20 increased a mean of 13 msec when load was applied (p < 0.05) 
and recovered following removal of the load. Mean APD 3 Q increased by 28 msec 
following the application of ventricular load (p < 0 . 0001 ) but only recovered partially 
by 1 min after removal of the load. Streptomycin. Streptomycin did not have a 
significant effect on APD and did not affect the response of APD 20 or APDg^ to 
changes in ventricular loading states. 



4 Model Analysis 

We investigated the possibility that altered intracellular calcium cycling via length- 
dependent activation may contribute to the increase in APD with stretch. Using a 
modification by Michailova and McCulloch (13) (Fig. 2) of the canine ionic model by 
Winslow et al (20), we simulated the effects of stretch by increasing the on-rate for 
calcium binding to troponin and studied the influence on action potential duration. 
This preliminary analysis showed that under some circumstances the resulting 
decrease in free myoplasmic calcium during systole was able to prolong relaxation 
somewhat (Fig. 3). 

Mechanoelectric feedback has been described in isolated cells and intact 
ventricular myocardium, but the mechanical stimulus that governs mechanosensitive 
channel activity in intact tissue is unknown. To investigate the effects of stretch- 
activated currents on regional action potential propagation in three-dimensions, we 
developed a finite element model of electromechanics in the rabbit heart. An 
anatomically detailed model of rabbit right and left ventricular geometry and fiber 
architecture was developed (15) and used to model regional mechanics during filling 
(16). See Fig. 4. The resulting strain field was then used as a stimulus for regional 
stretch-activated current (17). A stretch- dependent current was added in parallel to the 
ionic currents in the Beeler-Reuter ventricular action potential model (Fig. 5). 
Electrical propagation was simulated using the collocation-Galerkin finite element 
method. We investigated different mechanical coupling parameters to simulate 
stretch-dependent conductance modulated by either fiber strain, cross-fiber strain, or a 
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combination of the two. In response to pressure loading, the conductance model 
governed by fiber strain alone reproduced the epicardial decrease in action potential 
amplitude as observed in experimental preparations of the passively loaded rabbit 
heart (21). The model governed by only cross-fiber strain reproduced the transmural 
gradient in action potential amplitude as observed in isolated heart experiments, but 
failed to predict a sufficient decrease in amplitude at the epicardium. Only the model 
governed by both fiber and cross-fiber strain reproduced the epicardial and transmural 
changes in action potential amplitude similar to experimental observations (21). In 
addition, dispersion of action potential duration nearly doubled with the same model. 
These results suggest that changes in action potential characteristics may be due not 
only to length changes along the long axis direction of the myofiber, but also due to 
deformation in the plane transverse to the fiber axis. The model provides a framework 
for investigating how cellular biophysics affect the function of the intact ventricles. 




Fig. 2. Myocyte ionic model, adapted from (13) 




Fig. 3. Increasing troponin affinity for calcium prolonged the action potential in a preliminary 
model. Adapted from (12) 




Experimental and Computational Modeling of Cardiac Electromechanical Coupling 



117 




— = 0 S.cm^ 

= 5 MS.cm-5 




time (msec) 

Fig* 5* Action potential computed from Beeler-Reuter ionic model with varying degrees of 
stretch -activated conductance^ G^. Adapted from (1 1^ 17). 
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5 Discussion 

We have shown that acute loading of the left ventricle decreases apparent epicardial 
conduction velocity and increases action potential duration in the isolated rabbit heart 
and these responses do not appear to be mediated by a stretch-activated current. The 
decrease in conduction velocity could have potential arrhythmogenic consequences by 
decreasing the wavelength for reentry. In this way, mechanoelectric could contribute 
to arrhythmias associated with altered ventricular loading conditions. We also 
observed a prolongation of action potential duration with loading that was not blocked 
by streptomycin. A model analysis suggests that altered intracellular calcium cycling 
could be a contributing mechanism. Three-dimensional models provide a framework 
for integrating coupled cardiac functions across scales of organization from myocyte 
to whole organ. 



Acknowledgements. This research was supported in part by NSF grant BES-9634974 
and the National Biomedical Computation Resource, NIH P41 grant RR08605. 
Derrick Sung was supported by NIH training grant T32HL07444. Mary Ellen Thomas 
was an NSF scholar. 

Our modeling software can be downloaded free from http://cmrg.ucsd.edu/ 



References 

1. Dean JW, Lab MJ. Arrhythmia in heart failure: role of meehanieally indueed ehanges in 
eleetrophysiology. Laneet 1989;8650(1): 1309-13 12. 

2. Huang SK, Messer JV, Denes P. Signifieanee of ventrieular taehyeardia in idiopathie 
dilated eardiomyopathy: observations in 35 patients. Am J Cardiol 1983;51(3):507-512. 

3. Reiter MJ. Effeets of meehano-eleetrieal feedbaek: potential arrhythmogenie influenee in 
patients with eongestive heart failure. Cardiovaseular Researeh 1996;32(1):44-51. 

4. Hu H, Saehs F. Streteh-aetivated ion ehannels in the heart. J Mol Cell Cardiol 
1997;29(6):1511-23. 

5. Ruknudin A, Saehs F, Bustamante JO. Streteh-aetivated ion ehannels in tissue-eultured 
ehiek heart. Am J Physiol 1993;264:H960-H972. 

6. Craelius W, Chen V, el-Sherif N. Streteh aetivated ion ehannels in ventrieular myoeytes. 
Bioseienee Reports 1988;8(5):407-14. 

7. Sigurdson W, Morris C, Brezden B, Gardner D. Streteh aetivation of a ehannel in 
eollusean heart eells. Journal of Experimental Biology 1987;127:191-209. 

8. Salmon AH, Mays JE, Dalton GR, Jones JV, Eevi AJ. Effeet of streptomyein on wall- 
stress-indueed arrhythmias in the working rat heart. Cardiovase Res 1997;34(3):493-503. 

9. Eab MJ. Meehanoeleetrie feedbaek (transduetion) in heart: eoneepts and implieations. 
Cardiovase Res 1996;32(1):3-14. 

10. Sung D, Omens JH, MeCulloeh AD. Ventrieular meehanoeleetrie feedbaek in the isolated 
rabbit heart. In: Kamm RD, Sehmid-Sehonbein GW, Ateshian GA, Hefzy MS, editors. 
Summer Bioengineering Conferenee; 2001; Snowbird, Utah; 2001. p. 683-684. 

11. Riee JJ, Winslow RE, Dekanski J, MeVeigh E. Model studies of the role of meehano- 
sensitive eurrents in the generation of eardiae arrhythmias. J Theor Biol 1998;190(4):295- 
312. 




Experimental and Computational Modeling of Cardiae Eleetromeehanieal Coupling 119 



12. MeCulloeh A, Sung D, Thomas ME, Miehailova A. Chapter 11: Computational and 
Experimental Modeling of Ventrieular Eleetromeehanieal Interaetions. In: Virag N, Blane 
O, Kappenberger L, editors. Computer Simulation and Experimental Assessment of Cardiae 
Eleetrophysiology. Amonk, NY: Futura Publishing Company, Ine.; 2001. p. 89-94. 

13. Miehailova A, MeCulloeh AD. Modeling Ca^^ transients and Ca^^ and Mg^^ exehange with 
ATP and ADP during exeitation-eontraetion eoupling in ventrieular myoeytes. Biophys J 
2001;81(2):614-629. 

14. Sung D. Effeets of Meehanieal Load on Ventrieular aetion Potential Propagation and 
Repolarization [Ph.D.]. La Jolla: University of California San Diego; 2001. 

15. Vetter FJ, MeCulloeh AD. Three-dimensional analysis of regional eardiae funetion: a 
model of rabbit ventrieular anatomy. Progress in Biophysies and Moleeular Biology 
1998;69(2-3):157-183. 

16. Vetter FJ, MeCulloeh AD. Three-dimensional stress and strain in passive rabbit left 
ventriele: A model study. Ann Biomed Eng 2000;28:781-792. 

17. Vetter FJ, MeCulloeh AD. Meehanoeleetrie feedbaek in a model of the passively inflated 
left ventriele. Ann Biomed Eng 2001;29:414-426 (and eover). 

18. Sung D, Omens JH, MeCulloeh AD. Model-based analysis of optieally mapped epieardial 
aetivation patterns and eonduetion veloeity. Ann Biomed Eng 2000;28(9): 1085-1092. 

19. Sung D, Somayajula-Jagai J, Cosman P, MeCulloeh AD. Phase- shifting prior to spatial 
filtering enhanees optieal reeordings of eardiae aetion potential propagation. A nn Biomed 
Eng (in press) 2001. 

20. Winslow RE, Riee J, Jafri S, Marban E, ORourke B. Meehanisms of altered exeitation- 
eontraetion eoupling in eanine taehyeardia-indueed heart failure, II: model studies. Cire Res 
1999;84(5):571-86. 

21. Zabel M, Roller BS, Saehs F, Franz MR. Streteh-indueed voltage ehanges in the isolated 
beating heart: importanee of the timing of streteh and implieations for streteh-aetivated ion 
ehannels. Cardiovaseular Researeh 1996;32(l):120-30. 




Towards Model-Based Estimation of the Cardiac 
Electro-Mechanical Activity from ECG Signals 
and Ultrasound Images 



N. Ayache^, D. Chapelle^, F. Clement^, Y. Coudiere^, H. Delingette^, J.A. 
Desideri^, M. Sermesant*^, M. Sorine^, and J.M. Urquiza^ 



^ Epidaure Research Project, INRIA Sophia Antipolis, 
2004 route des Lucioles, 06902 Sophia Antipolis, France 
^ Sinus Research Project, INRIA Sophia Antipolis 
^ Macs Research Project, INRIA Rocquencourt 
^ SOSSO Research Project, INRIA Rocquencourt 



Abstract. We present a 3D numerical representation of the heart which 
couples the electrical and biomechanical models. To achieve this, the 
FitzHugh-Nagumo equations are solved along with a constitutive law 
based on the Hill-Maxwell rheological law. Ultimately, the parameters 
of this generic model will be adjusted by comparing the actual patient’s 
ECG with computational results and the deformation of the biomechani- 
cal model with the geometric information extracted from the ultrasound 
images of the patient’s heart. 



1 Introduction 

The objective of our multidisciplinary project ICEMA (standing for Images of 
the Cardiac Electro-Mechanical Activity, a collaborative research action between 
different INRIA projects) is to build a generic dynamic model of the beating 
heart and a procedure to automatically adjust the parameters to any specific 
patient. We plan to construct an identification procedure using 2 sets of rela- 
tively easy-to-access measurements on a patient: the ECG (Electrocardiogram), 
and a time sequence of volumetric ultrasound images. Once the generic model 
is adapted to a specific patient, it becomes possible to derive a set of quan- 
titative and objective parameters useful in helping clinicians and physiologists 
to better understand the electro-mechanical coupling and diagnose pathological 
conditions. Significant results are expected in the following fields of cardiovas- 
cular pathology: 

— assessment of the haemodynamic repercussions of heart rate and electrical 
conduction disorders; 

— assessment of the degree of heart failure (inefficiency of the cardiac pump 
induced by weakened cardiac contractility resulting in increased local wall 
thickness and decreased ejection fraction); 
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— assessment of the electrical and mechanical repercussions of cardiac infarc- 
tion (right ventricle infarction is believed to rather lead to electrical reper- 
cussions while left ventricle infarction is believed to lead to mechanical reper- 
cussions); 

Our approach combines a 3D numerical model of the electric wave prop- 
agation with a 3D biomechanical model of the cardiac muscle. The 2 models 
are explicitely coupled in the simulation to generate a dynamic behaviour of 
the heart. The model for electric wave propagation is derived from FitzHugh- 
Nagumo equations, while the mechanical model is based on the classical Hill- 
Maxwell rheological law. These models are expected to reflect on a macroscopic 
scale the coupling present on the cellular scale. To provide a realistic motion of 
a standard beating heart, the highly anisotropic nature of the muscle fibres in 
standard anatomy are accounted for. The electric wave is propagated from the 
extremities of the Purkinje network. 

Two error functions will serve to adjust the parameters of this generic model 
to a specific patient: the first will compare the actual patient’s ECGs with a set 
of ECGs computed from the simulation. The second will compare the deforma- 
tion of the biomechanical model with the motion extracted from the ultrasound 
images of the patient’s heart. In addition, the ventricular blood pressure, when 
available, is readily introduced as a constraint for the deformation of the biome- 
chanical model. Ultimately, a retro-action procedure will be used to update the 
parameters of the generic model from these error functions. 




In this article, we develop the current stages of this on-going research. In 
section 2, we describe the construction of a finite element model incorporating 
standard anatomical knowledge. In sections 3 and 4 we respectively describe the 
electrical and mechanical models. Section 5 describes a preliminary approach 
to account for the geometric information provided by the ultrasound image se- 
quence. Finally, section 6 discusses the possible refinements we plan to introduce 
in each of these models, and the directions towards the identification of the pa- 
rameters of the global system. 
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2 Anatomical Model 

2.1 Geometry and Muscle Fibres 

Our geometric model consists of a mesh of the ventricular myocardium (including 
both right and left ventricles) and of the myocardium fibre directions. Indeed, 
the fibre architecture has a great influence on the motion and propagation of 
the electrical excitation. The model is based on data available from the Bioengi- 
neering Research Group of the University of Auckland, New Zealand [9]. These 
data have been obtained by the dissection of a dog’s heart and are composed of 
a mesh with 256 nodes, 180 hexahedra, and fibre direction at each node. 

Since these data are too coarse for 3D computations, new meshes have been 
generated using a procedure which consists in: 

1. interpolating the surface of the mesh with a mesh of triangles; 

2. remeshing the 3D volume with tetrahedra; 

3. interpolating the fibre directions on the new mesh. 

We obtained meshes with 623, 3763, and 9623 vertices. 

2.2 Conduction Network and Electrical Onset 

A complete three-dimensional model of the heart should include both the car- 
diac muscle and special conduction network (Purkinje fibres and His bundle 
branches) [2,13] and whose role is rendered through boundary conditions: the 
fibres of Purkinje are assumed to end up within the muscle in a region of the en- 
docardium; a pulse-shaped boundary condition is applied on this hand-delimited 
region as a model of onset of the cardiac excitation. 



3 Electrical Model 



Each heart-beat cycle, a depolarisation wave is initiated (see above) and propa- 
gates along the myocardium during about 10% of the total cardiac cycle. 

Among the various models for the electric wave propagation, the FitzHugh- 
Nagumo model is classical [1]. It writes: 



— = div(L>V'u) + f{u) - V, 

dv , . 

— = e{ku - v). 



( 1 ) 



where u represents the transmembrane potential and v is an auxiliary variable, 
and f{u) = fQu{l — u){u — a). The system is subject to the boundary condition: 

u{x,t) = e{x,t) X e 



where e is the pulse-shaped input signal. 
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Anisotropic behaviour is accounted for by making the diffusion matrix D 
depend on the fibre orientation: Vx, D{x) = diag{d^s^s) in the orthonormal 
basis (n, k, 1) local to x ; where n is a unit vector tangent to the fibre. 

The theoretical aspects of (1) have been widely studied [16]: a travelling 
wave of fixed shape and speed either appears or not, depending on the initial 
excitation being above or below a threshold. The important parameters are 

— a which controls the value of the threshold for excitation; 

— D and /o which adjust the wave speed and time scale. 

A simple output, the electrocardiogram along the axis d, can be calculated 
as follows: 



ecg(d, t) = L{\/u{t), d). 

1/ is a weighted integral over the volume of the heart. 

Although more accurate models exist [1], we retained this one, as it correctly 
captures in a simple formulation, the qualitative behaviour of the wave and leads 
to fast three dimensional computations [14,10,7,13]. This model would likely be 
improved by introducing an additional variable [8]. 

Solutions to (1) have been approximated by using a standard PI Lagrange fi- 
nite element procedure (with mass lumping and first order numerical integration 
at vertices), on the given tetrahedral anatomical mesh. Since the time-dependent 
phenomenon is of interest, small time-steps are used in an explicit procedure 
(Euler method). 

The parameters a, e, /c, /o are estimated according to the results in [8]. 
Potential maps such as those illustrated on Figure 1 are derived at different 
time steps, starting after a wave has been initiated at the apex and propagated. 
Better initiation points will be chosen according to data from [5]. 

Fig. 1. (Top) Isotropic propagation (d = 1, e = 1), (Bottom) Anisotropic propagation 
(d= 1, £ = 0.7) 
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This time-dependent computed potential can then be used as an excitation 
entry to the system describing the mechanical behaviour of the myocardium. 
No particular feedback is currently used, as proposed in [8] to strengthen the 
coupling between the electrical and mechanical models. 

4 Mechanical Model 

The constitutive law that we consider for the cardiac muscle (composed of a 
distribution of stacked myocardial fibres) is based on the classical Hill-Maxwell 
rheological model [6], schematically represented in Figure 2 where Ec^ Eg and 
Ep symbolically denote the contractile, series and parallel element, respectively. 




Fig. 2. Hill-Maxwell rheological model 



The series and parallel elements are elastic, but non- necessarily linear. The 
contractile element develops an action oriented in the direction of the fibre at 
the point in consideration, hence the corresponding stress tensor is everywhere 
one-dimensional (ID) and in the form ac = (Jdk^E^ where n denotes the unit 
vector tangent to the fibre. This implies that the stress tensor corresponding to 
the series element is also ID, since = ^c- In the parallel element, the stress 
tensor is fully 3D (even though the corresponding constitutive law is strongly 
anisotropic due to the fibre architecture), so is the global stress tensor given by 
E = Ec + Ep- 

For the (scalar) stress quantity cFc we use a (time-dependent) constitutive 
relation recently proposed and justified in [3] (see also [4]), viz: 

dc = - (|4| + d\u\) dc + kcSc + CTo|i4| + , 
kc = — (|^c| + d\u\) kc + /co|'u| + , 

(Jc = kcSo -h dc + 1^6 c. 

Here, 5c denotes the (ID) strain corresponding to the contractile element, kc 
and dc represent additional internal variables, u is the input (also represented in 
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Fig. 2), and all other quantities are constants. The input is a chemical quantity 
that depends mainly on the calcium concentration [4] , and which is modelled in 
section 3. Hence, the resulting constitutive law for the global rheological model 
is of the form = U {e, u), where U denotes a functional. This constitutive 
law is then used in combination with the equations of solid dynamics to obtain 
the mechanical model. 

5 Interaction with Ultrasound Images 

In our current implementation, we simplify the mechanical model into 
anisotropic linear elasticity and the fibres activation is modelled as a stress ten- 
sor a a = / (8) /, where a is the activation rate which is directly related to the 

potential u of Section 3, and / the fibre direction. It gives a force F^: 

Fa= div{aa) dv = / div{af 0 f) dv = / {af0f)0nds. 

Jv Jv Js 

Therefore, when a fibre is activated, its contraction is equivalent to a pressure 
applied to the surface of the tetrahedron in the fibre direction. 

Figure 3 presents the results of a contraction proportional to the potential 
(which should be replaced by a differential equation controlling the activation 
from the potential). 




Fig. 3. Effect of the fibres contraction on the 3D model. 



The model is also constrained by the ultrasound images (see Fig. 4) through 
the external forces F^ which are proportional to the distance to the closest bound- 
ary point of the image from the considered point of the mesh. 

We use mass-lumping in a Newtonian differential equation with an explicit 
integration scheme to compute the position P of each vertex: 



2At J 



>t+i _ 



= Fi 



— P* - ( — -A- 



At^ 



2 At) 



with 7 the damping factor. F represents the forces computed from linear elas- 
ticity plus the boundary conditions (in particular the ventricular pressures). As 
ultrasound data is sparse and noisy, it is efficient to use a volumetric model [11, 
12] together with electro- mechanical knowledge to extract information (see [15] 
for preliminary results). 
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Fig. 4. Left: the 3D model in a slice of the 3D ultrasound image. Right: intersection 
of the model and the image. 



6 Perspectives 

Each part of our model is planned to be improved: 

— the anatomical model (Section 2) should be based on a human heart. Diffu- 
sion MRI is probably one possible way to obtain fibre directions and maybe 
Purkinje network; 

— the electrical model (Section 3) could include a third variable to better con- 
trol the shape of the wave and a mechano-electrical feedback; 

— the mechanical model (Section 4) has to be solved on a 3D mesh; 

— the simplified model for ultrasound segmentation (Section 5) will become 
non-linear. 

For each of these models we intend to identify parameters: 

— for the electrical parameters, D, /o, and the electrical entries, Purkinje fi- 
bres, ectopic foci, should be estimated by comparing computed ECGs with 
measured ones. To achieve this, an inverse problem has to be considered; 

— for the mechanical parameters, identification techniques still have to be de- 
veloped, as in vivo rheological studies for human tissues are hard to set up. 
Here, a criterion will be the difference between the computed motion and 
the one extracted from the ultrasound images. 

All these points will be the topics of our future work. 
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Abstract. We present a 3D mechanical model for the behaviour of the 
cardiac muscle, based on a recently proposed model of myofibre con- 
traction. The construction of the 3D model involves a rheological model 
similar to that of Hill-Maxwell. We then introduce some discretisation 
techniques adapted to our mechanical model, and we report on some 
preliminary numerical results. 



1 Introduction 

Like most biological materials, the cardiac muscle displays strongly nonlinear 
and anisotropic mechanical behaviour, see, e.g., [5]. Anisotropy is governed — in 
particular — by the architecture of muscle fibres and of the fibre sheets. Moreover, 
muscle fibres are responsible for the active stress generated by the electrical 
stimulations. This active stress mainly concerns the stress component along the 
fibre direction, which is also the primary direction of electrical propagation. 

In this work we present a complete mechanical model that aims at cor- 
rectly representing the cardiac muscle behaviour in the presence of the chemi- 
cal/electrical stimulation. This model is based on a recently proposed myofibre 
model [2], which we recall before presenting the methodology allowing its in- 
corporation into a global three-dimensional (3D) mechanical model. The last 
section is devoted to numerical discretisation strategies adapted to the model, 
and we present some preliminary numerical results obtained for a ID simplified 
problem. 

2 The Excitation-Contraction Myofibre Model 

The myofibre model introduced and justified in [2] (see also [3]), is based on 
the classical Hill-Maxwell rheological model [6], depicted in Figure 1. Elastic 
material laws are used for the series element Eg and for the parallel element E^. 
Based on empirical results, the corresponding stress-strain laws are generally 
assumed to be of exponential type [9,12]. 

The element Ec accounts for the contractile electrically-activated part of 
the behaviour. In order to represent this behaviour, we use the constitutive law 
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recently proposed in [2] (a refinement of the model presented in [3]), based on the 
microstrnctnral sliding filament model of Huxley [8] and the distribution-moment 
approach of Zahalak [13]. This (time-dependent) constitutive law reads 

{ CTc = -(k'cl + \u\),ac + kJc + (Jo\u\+, 

kc = -{\ec\ + \u\)kc + ko\u\+ ( 1 ) 

(Tc = kcEo + CTc + I^Sc, 

where u is the excitation, a chemical quantity that mainly depends on the calcium 
concentration, 5c stands for the time derivative of 5c, and ao, ko, 5q, n are local 
parameters such that ao, /cq, > 0 . This model is also in accordance with 
observations made on the passive behaviour of the muscle [10]. 




Fig. 1. Hill-Maxwell rheological model 



3 Geometry, Fibre Orientations, and Governing 
Equations of Activated Heart Contractions 

In order to construct our global 3D mechanical model we need: 

— Anatomical data allowing for a definition of the heart geometry, with an 
accurate description of the microstrnctnral anatomy, i.e. of the orientations 
of muscle fibres and fibre sheets. 

— The governing equations of the 3D mechanical behaviour. 

We used the anatomical data provided by the Auckland biomedical engineer- 
ing research group [7]. We refined the original mesh composed of hexahedral 
elements by using the softwares YAMS and GHS3D developped at INRIA in 
the Gamma team [4,11]. The refined mesh, composed of tetrahedral elements, is 
shown in Figure 2. 

The orientation of the fibres is available at each integration point in the mesh. 
The corresponding streamline representation is shown in Figure 3. 
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Fig. 2. Refined tetrahedral mesh 



Fig. 3. Heart fibre directions 



We now introduce the governing 3D mechanical equations. We denote by 
F and s the deformation gradient and the Green-Lagrange strain tensor, re- 
spectively. We then use the second Piola-Kirchhoff stress tensor, denoted by a. 
Recalling the rheological model of Figure 1, the total stress is the sum of the 
stresses in the two branches. The branch containing the contractile element, 
however, only contributes to the stress tensor in the direction of the fibre, given 
by the unit tangent vector n. Namely, we have 

a = ap + (JiD G F, (2) 

where (Jid is a scalar quantity that represents the (second Piola-Kirchhoff) stress 
in the series branch. The series kinematics and the one-dimensional character 
of the deformation in this branch imply that the corresponding (scalar) Green- 
Lagrange strain 

SiB = '^eijniJij (3) 

ij 

satisfies the multiplicative law 

1 + ^ID = (1 + ^c)(l + ^s): (4) 

where 5c and Ss denote the strains in the contractile and series element, re- 
spectively (note that this relation reduces to the usual additive law in the case 
of small strains). Furthermore, the associated second Piola-Kirchhoff stresses 
satisfy, by formal thermodynamical considerations, 

1 Ss 1 Sc 



= 



( 5 ) 
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On the other hand, the parallel kinematics directly implies 

g = (6) 

Combining the above equations with the constitutive laws of the elastic elements 

Gs — Gg(Sg), (7) 

and with the equations (1), the 3D material behaviour is completely defined. 

The global 3D mechanical model is obtained by using these constitutive re- 
lations (involving the input u) together with the equation of dynamics 



div(£ ’g) - py = 0, (8) 

where y denotes the displacement vector and p the mass density (in the original 
configuration) . 



4 A One-Dimensional Activated Fibre Contraction 



As a model problem we consider a one-dimensional homogeneous medium gov- 
erned by a ID reduction of equations (l)-(8). 

In this first step to the numerical exploration of the excitation- contract ion 
model we use the infinitesimal displacement assumption: with y{x,t) denoting 
the longitudinal displacement of a material point located at position x, 0 < x < 
1, in a reference stress-free state, the deformation is given by e{x^t) = 
Moreover, we assume that series and parallel elements are linear and that v = 
0, 5o = 0. The considered equations read 

' Py- + = 0 , 

be = -(141 + |^|)cTc + A:c4 + cro|^| + , /qn 

kc = -(| 4 | + \u\)kc + ko\u\^, 

^ kg(^£ 



with positive constants kp and kg. We also consider that the body is fixed at 
both ends, that is ^(0, t) = y{l,t) = 0, corresponding to an isometric contraction 
or relaxation. 

Performing a variational finite element space discretization with piecewise- 
linear continuous functions for the displacement and piecewise-constant func- 
tions for stress and deformation variables, we use a midpoint scheme for solving 
the resulting system of differential-algebraic equations (see [1]). 

In the numerical tests presented in Figures 4 to 7 we chose all initial un- 
knowns to be zero. These numerical experiments illustrate the active contraction 
(positive u) and active relaxation (negative u) abilities of the model. 
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Fig. 4. Electric excitation 



Fig. 5. Resulting stress 



yixo 





Fig. 6. Deformation 



Fig. 7. Displacement 
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Abstract. We present a one-dimensional mathematical model to ex- 
plain changes in regional radial myocardial deformation during is- 
chaemia. The model makes use of cellular contraction properties of nor- 
mal and ischaemic myocytes extracted from an animal model of induced 
chronic, regional ischaemia of the posterior wall. Simulations of the model 
are compared to the deformation patterns observed in the experimental 
animal models of regional ischaemia. 



1 Introduction 

Localized, regional ischaemia induces reduced systolic thickening with a concomi- 
tant increase in post-systolic thickening (PST). Prior experimental and clinical 
observations proposed acute development and presence of PST in an ischaemic 
myocardial segment as a potential marker of myocardial viability. With the in- 
troduction of (echocardiographic) strain and strain rate imaging as a new clinical 
tool to non-invasively quantify regional myocardial deformation [1], evaluating 
explicit markers such as PST can increase the clinical potential of the method. 
Controlled experimental work, e.g. [2, 3, 4, 5], has shown a complex mechanical 
behaviour of the ischaemic myocardium. In summary in ischaemic myocardium 
the following observations can be made: (a) during the isovolumetric contrac- 
tion period and early systole, myocardial thinning occurs, (b) in systole, the 
magnitude of regional thickening/shortening is decreased, (c) after aortic valve 
closure, normal relaxation and thinning (lengthening) is interchanged with on- 
going post-systolic thickening (shortening). The magnitude of these alterations 
is proportional to the severity of the ischaemic insult. This phenomenon occurs 
in acute, prolonged as well as in chronic regional ischaemia. Despite the multi- 
tude of experimental work, the mechanisms underlying PST remain unclear. The 
debate is still going on whether this phenomenon is active, i.e. the thickening 
in the ischaemic segment is due to late, prolonged contraction in the ischaemic 
fibres, or passive, thickening due to the elastic force from the interaction be- 
tween the relaxing normal segments and the ischaemic segment. In this paper. 
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we try to make a further step in the explanation of this phenomenon. Assum- 
ing that PST mostly depends on regional heterogeneity in contractile and/or 
elastic properties of interacting segments of myocardium, we have developed a 
simple one-dimensional model reflecting these properties, which explains these 
phenomena, at least qualitatively. 

2 Theoretical Model 

In this section we present a simple model of normal and ischaemic myocardium, 
interacting during a cardiac cycle. 



2.1 Description 

Though, anatomically accurate models have been constructed for the heart, we 
will rely on a simplified model, using finite difference techniques to solve the 
dynamics. Despite the fact that this model has severe limitations, it is able to 
simulate qualitatively the regional deformation patterns of the myocardium. The 
choice of the simple model enables us to solve the equations of motion for the 
nodes in fig. 1 straightforwardly, by simple time step integration techniques, in 
a limited computation time (±50s), using Matlab (The MathWorks, Inc.). 
Model geometry. In the model a midwall segment of the left ventricle is effec- 
tively described by a discretised, closed chain where the thickness of the wall is 
defined as a field on this chain. The reduction from the three dimensional struc- 
ture to the one-dimensional model is depicted in Fig. 1. This effective model takes 
out the interactions in the longitudinal directions of the left- ventricular wall. The 




Fig. 1. Representation of the reduction steps from a mid- wall slice of myocardium to 
the discrete one-dimensional description (red: ischaemic, blueinormal) 



model allows for a connected variable sized region of ischeamic segments along 
the chain. The “thickness” of the elements is defined to be proportional to the 
reciprocal of the length of the elements (based on local mass conservation). In 
fig. 1, half of the chain is ischaemic and half is normal. The elements are described 
as Maxwell elements splitting active and passive properties of the elements. It 
describes the element as a contractile component (CE) in series with an elastic 
component (SE), together in parallel with another elastic component (PE). 
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Passive elastic properties. The SE and PE are described as non-linear springs, 
with a passive length-tension relation, i.e. 

T = - 1) , where A = , (1) 

l{t) is the instantaneous length and Iq the end-diastolic length of the segments. 
The parallel component of the Maxwell elements is assumed to take into account 
the elasticity of the extracellular collagen matrix. 

Active properties. The active developed contraction force of the segments is 
modeled by a gaussian, with variable amplitude and width, which can be shifted 
in time, allowing variation of the amplitude, duration and onset of the active 
contraction. Again, this is an extensive simplification of the active contraction 
sequence. Complete models includes a large number of parameters and state 
variables, which leads to a better representation but also decrease the possibility 
of estimating parameters from them. 

Boundary Conditions. As a boundary condition the intercavitary pressure was 
simulated according to measured pressure curves (fig. 3). 



2.2 Simulations 

To simulate the radial deformation during a complete cardiac cycle, the cycle 
was divided into 800 time steps of 1 ms. We divided the cycle into 4 events, 
isovolumetric contraction period (0-50ms) (end: aortic valve opening), ejection 
(50-300ms) (end: aortic valve closure), isovolumetric relaxation period (300- 
350ms) (end: mitral valve closure) and diastole (350-800ms) (end: mitral valve 
opening). Equations of motion were integrated at each time step. The chain was 
divided into 50 segments /nodes and had a diastolic radius of 2 cm, corresponding 
to an average diastolic cavity radius of the pig model of section 3. 

During the isovolumetric periods motion of the nodes perpendicular to the 
instantaneous segment directions was restricted. This perpendicular motion was 
not restricted to zero because deformation and motion caused by shape changes 
can occur during the isovolumetric periods. 

The final result of the simulations was the natural transmural radial strain, 
which represents wall-thickening during the heart cycle. 

We have performed three groups of simulations to study myocardial is- 
chaemia: reduced perfusion, total occlusion and infarction. Eor reduced perfusion 
and total occlusion, the passive properties of the ischemic zone were taken to 
be normal, but in the case of reduced perfusion, the active contraction force 
was prolonged and reduced, in accordance with the results of cellular contrac- 
tion measurements from the experimental pig model, (see section 3). Eor total 
occlusion the active force equals zero. Infarction was simulated by zero active 
force, but with higher stiffness constants in the ischaemic zone. The qualitative 
simulation results for radial strain are summarized in fig. 2. 
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(A) 





(B) (C) 



Fig. 2. Simulation results. (A)Reduced perfusion. (B)Total occlusion. (C)Infarction. 
Dashed line: Radial strain in the middel segment of the normal zone, Solid line: Radial 
strain in the middle segment of the ischaemic zone. AVC: Aortic valve closure, MVO: 
Mitral valve opening. 



3 Comparison of the Model with Experimental Data 

In this section we shortly describe a pig model for induced regional chronic 
ischaemia and present some data. This closed chest model was developed to 
get a spectrum of chronic regional ischaemia, reflecting the clinical setting. The 





Fig. 3. Experimental data. (A) Transmural radial strain (wall thickening) curves from 
the chronic ischaemic posterior wall of the pig model, (B) Transmural radial strain curve 
from a pig experiment of acute total occlusion [5]. (C)Pressure curve measured in exper- 
iment and used as boundary condition for simulations. (D) Single twitches from stim- 
ulated single cells (dashed: normal myocytes, solid: ischaemic myocytes). AVO/AVC: 
Aortic valve opening/closure, MVO/MVC: Mitral valve opening/closure. 



model has been developed, by implanting a copper coated stent in the proximal 
segment of the left circumflex coronary artery. Two weeks after implantation 
the percent diameter of the developed stenosis was estimated by angiography 
and ultrasound strain rate imaging was performed according to the protocol 
described in [5] . After sacrifice of the pigs the heart was taken out and prepared 
for cell extraction. The radial strain curves of the posterior wall are depicted 
in fig. 3(A). Cells of the ischaemic region of the excised heart were isolated 
and compared with cells of the same region of normal hearts. Cell contraction in 
isotonic twitches of single cells have been measured and are depicted in Fig. 3(D). 
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4 Discussion 

The proposed simple one-dimensional model is able to simulate the regional 
radial deformation curves for a whole spectrum of ischaemic myocardium, as 
observed in diverse experimental settings [4,5]. The model is based on only a 
small number of parameters for normal and ischaemic regions, i.e. elasticity 
constants together with strength, length and onset of actively developed force. 
From the close resemblance between measured and simulated deformations we 
can infer that these are the major determinants of the observed phenomena. 

One of the major objectives of the simulation study was to get insight into 
the physiological mechanisms behind PST, an abberant myocardial thickening, 
occuring in what is traditionally perceived as the relaxation (and thus thin- 
ing) phase of the cardiac cycle. Two different explanations have been proposed, 
namely “active” and “passive” thickening. Active thickening after aortic valve 
closure could occur if there is a delayed and/or prolonged active contraction 
force developed in the myocytes. On the other hand, passive thickening could 
be explained by the interaction with neighbouring segments. In this hypothesis, 
ischaemic tissue would not thicken by itself, but exhibits a pulling force from 
the neighbouring thickening (normal) myocardium. Since, before aortic closure, 
pressure is high in the cavity, this pulling force is not sufficient to deform the 
non-contracting part as much as the normal tissue, resulting in less deformation 
(explaining the observed decrease in systolic thickening). However, after the sud- 
den drop in pressure associated with aortic valve closure, the pulling force from 
the normal (and at that moment thicker) tissue can further deform the ischaemic 
region, until the dimensions are similar to the (at that moment thining) normal 
neighbouring part. This would explain the observed PST. 

The implemented model was developed in such a way that both hypothe- 
ses could be investigated. Therefore, elastic properties and myocyte contraction 
parameters (for the actively developed forces) were combined. 

To investigate whether the passive theory could completely explain PST in 
all cases, a total occlusion simulation, i.e. no active force is present, was per- 
formed (Fig. 2(B)). When comparing findings with experimental total occlusion 
(Fig. 3(B)), the observed deformation patterns are very similar, confirming the 
passive explanation. Also, when tissue elasticity is decreased in the model , mim- 
icking the development of fibrosis after infarction, myocardial strain could be 
predicted by the simulation (Fig. 2(C)), solely based on passive phenomena. 

However, strain curves, measured during experimental reduced perfusion, 
cannot be explained by the model when only incorporating passive deformation. 
From the experiments in isolated myocytes, we could observe that in chronic 
ischaemia there is a prolonged, but reduced contraction of the cells, with no 
delay in onset. Consequently, this observation was used as input for further sim- 
ulations. Fig. 2(A) shows the result of the simulation where, additionally to 
the passive interaction of segments, a reduced and prolonged contraction force 
was imposed on the ischaemic region. Under these conditions, it can be clearly 
observed when comparing to the measured strain curves in case of experimen- 
tally reduced perfusion (Fig. 3(A)), that this combination of the passive and 
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active hypotheses explains both the systolic and the post-systolic deformation 
of ischaemic myocardium. 

Additionally, in all cases, the simulations exhibited early thinning of the 
ischaemic segments. During the isovolumetric contraction period the segments 
can rearrange in the circumferential directions. During this early period in the 
cardiac cycle, the normal segments develop a higher contraction force and stretch 
the ischaemic segments, which results in a temporary thinning. 

As a conclusion, we can state that our model would support the following 
hypothesis regarding myocardial deformation during ischaemia: With decreasing 
coronary flow, there is a decrease in systolic thickening caused by a reduction of 
the active contraction force developed in the myocytes. After aortic valve closure, 
there is a further thickening (PST) where a contineous thickening is caused by 
a delayed active contraction and sudden changes in thickness can be explained 
by passive interaction with surrounding tissue, after all active force development 
ceased. The more the perfusion is diminished, the more the deformation will be 
governed by passive mechanisms. 

Finally, we remark that although the model explains qualitatively the mech- 
anisms of PST, it has its limitations. It does not take into account longitudinal 
deformation. Additionally, in the diastolic phase the model can not explain the 
biphasic responses of the radial strain. In part this is due to the fact that we 
take the overall pressure as the external force that act on the myocardium, while 
in reality this strain is influenced by the propagation of the pressure wave and 
its associated vortices. Therefore, quantitative results from the model cannot be 
trusted. 
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Abstract. We deseribe simulations of propagated eleetrieal exeitation in three- 
dimensional anisotropie myoeardial musele. Aeeording to the bidomain theory, 
anisotropie eleetrieal eonduetivities are presented as tensors in the intraeellular 
and interstitial domains {Di and De, respeetively). Under the assumption of equal 
anisotropy ratio (Di = kDe), subthreshold behaviour of the exeitable elements is 
governed by a parabolie reaetion-diffusion equation for the membrane potential, 
solvable even on a desktop eomputer. In the ease of more general anisotropies 
(Di / kDe), also the interstitial potential needs to be solved simultaneously 
from an elliptie partial differential equation, requiring a supereomputer for large 
arrays of exeitable elements. In both eases, the elements obey eellular automata 
rules in the suprathreshold state. We present preliminary results of the propagated 
exeitation for different anisotropy ratios in a three-dimensional slab geometry. 



1 Introduction 

First idealized models describing the normal activation sequence in the human heart 
were reported over two decades ago [1,2,3]. Mainly due to computational limitations, 
the models did not include myocardial anisotropy nor physiological propagation. On 
the basis of anisotropic bidomain theory and cellular automata theory, development of 
more realistic whole-heart models have become feasible. A ventricular model that pro- 
duces a correct normal activation sequence is a prerequisite for simulating pathological 
conditions, such as ischemia, infarction or ventricular arrhythmias. 

We describe here the theoretical basis of the Dalhousie hybrid propagation model 
[4,5,6], based on equal anisotropy ratio and cellular automata. We present how the 
model can be extended to more general anisotropies and show our first simulations in 
a three-dimensional slab geometry. Due to extensive computational requirements, such 
simulations may be too demanding to cover the whole ventricles 2 x 10^ elements), 
but they are still useful to explore the validity of the simpler equal anisotropy approach. 
Furthermore, unequal anisotropy may be necessary to simulate detailed characteristics 
of the epi- and endocardial potential distributions. 
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2 Bidomain Theory 

The basis of our model is the bidomain theory for anisotropic myocardium, which rests 
on the assumption that the cardiac tissue functions as an electrically conductive syn- 
cytium consisting of two interpenetrating domains — intracellular (i) and interstitial (e) — 
connected everywhere via the cell membrane [2,3]. Consider the anisotropic bidomain 
region H, which represents a system of cylindrical myocardial fibers with anisotropic 
electrical conductivity [7,8]; H is embedded in a bounded, insulated and isotropic con- 
ducting medium B.A tissue structure in H is assumed to have longitudinal fiber direction 
defined by a unit vector (x) that is allowed to vary with position x, while intracellular 
and interstitial domains share the same (x) . The local basis consisting of an orthogo- 
nal set {ai(x), a2(x), a3(x)} is chosen at x so that a unit vector a3(x) is parallel with 
a^ (x) ; unit vectors ai (x) and a2 (x), and all vectors coplanar with them, are said to be in 
the transverse direction at x. In the local basis, the intracellular conductivities along the 
axes are cr|, al, erg, and the corresponding interstitial conductivities are af , af, cr|; due 
to axial symmetry of each fiber, ^ are conductivities in the transverse 

direction, and dg ^ are conductivities in the longitudinal direction. Moreover, 

uniform anisotropy assumption requires that scalar constants are independent 

of X. 

2.1 Conductivity Tensors 

Let us introduce second-rank tensors D* and D* to describe comprehensively the intra- 
cellular and interstitial anisotropic conductivities [8]. In the local basis the conductivity 
tensors are diagonal: D* ^ = diag(cr^’^, To allow variable fiber direction, 

set the global (cartesian) coordinate system in which local basis is defined at any x as 
A = (ai(x), a 2 (x), a 3 (x)), where ai(x), a 2 (x), and a 3 (x) are the column vectors. To 
represent tensors B* and D* in the global coordinate system, the local basis has to be 
rotated by multiplying it on each side by a rotation matrix. In the global coordinate sys- 
tem the conductivity tensors are then where the superscript T denotes 

matrix transpose. Therefore, we can write the conductivity tensors Bi^e as [6] 

A,e = (crf-aj->,af + aj-V, (1) 

where / is the identity matrix. Thus, according to Eq. I, an anisotropic medium can 
be thought of as having isotropic conductivities throughout, plus the “boost” in 
conductivities, (a^’^ — along the fiber direction. In addition to conductivity tensors 
Bi and Be , we will introduce a conductivity tensor B = BiABe^ which will characterize 
the composite medium of the bulk cardiac tissue. 

2.2 Electrical Potential Distribution and Current Flow 

Electrical potential and current density are defined in H as macroscopic quantities, which 
may be considered to be averages over small volumes of cardiac tissue encompassing 
several cells. The current densities are given by the Ohm’s law: 

Ji = -BiV(j)i, Je = -BeV(t)e in H, 

Jo = -c^oV^o in 



( 2 ) 
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where and 0e are, respectively, electrical potentials in intracellular and interstitial 
domains; 0o is the extracardiac potential and (Tq is the isotropic scalar conductivity 
outside of H. In the bidomain region, the total macroscopic current density, J, is then 

J = Ji + Je = - DeV^e in (3) 

It is well established that the electrocardiographic volume-conductor problem can 
be treated as a quasi-static one [9]. Under the quasi-static assumption, the conservation 
law requires that the divergence of total current density vanishes (V • J = 0), which 
leads to 

V • ( A + A V0e) =0 in (4) 

i.e., current can flow from one domain to the other, but there are no net sources or sinks 
of current in H. Moreover, since thereare no sources in B, 

V • Jo = V • ctqVA = 0 in 5. (5) 

The intracellular and interstitial domains are coupled through a distributed cellular 
membrane; the outward transmembrane current per unit volume is, in accordance with 
the conservation law (Eq. 4), 

A = V • AV0i = -V • AV0e • (6) 

The transmembrane potential is deflned in H as 

~ • C7) 

The terms DiV(j)e and —DiV(j)e can be added to Eq. 3 and, by using Eq. 7 and recalling 
that D = Di ^ get for the total current density 

J = -DiVVm - DV(t)e = J' - DV(t)e • ( 8 ) 

Here J* was substituted for —DiVvm\ this term has a dimension of a current dipole 
moment per unit volume and it is called an impressed current density [9]. The impressed 
current density is driven by the electrochemical generators in cardiac cells, while the 
term —DV(j)e represents passive return currents in the tissue. From the condition that 
the divergence of J must vanish, or alternatively, by substitutions from Eq. 7 into Eq. 6, 
we can obtain a fundamental partial differential equation in (j)^, with a source term that 
involves the gradient of the transmembrane potential (c.f. Geselowitz and Miller [7] 

Eqs. (4) and (7); Plonsey and Barr [10] Eq. (7)) 

V • DV(t)e = -V • D,Vvm = V • JU (9) 

Thus, the distribution of the interstitial potential, (j)e, is related by Eq. 9 to the current 
sources. Eq. 9 suggests that H can be regarded as a composite medium characterized by 
the bulk conductivity tensor D, in which there is a distributed impressed current density 
(a current dipole moment per unit volume), deflned by the term J* = —DiVvm', the 
divergence of the latter is the impressed scalar current per unit volume [9]. 
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From the cable theory [1 1] we have 









dVr, 



dt 



lion la 



(10) 



where x is the membrane surface area per unit volume of the tissue, Cm is the membrane 
capacitance per unit area of the membrane surface, Con is the ionic current per unit area 
of the membrane surface, and lapp represents an applied current stimulus to start the 
activation. By combining Eqs. 6 and 10 and by defining 



xa 



X^ic 



^app 



xi. 



app 5 



we get an equation that relates the spatial distribution of intracellular potential and the 
membrane dynamics: 

dVm 

Cm ^ Hon ^app ~ ^ V (/)i iu H . (11) 

By substituting = Vm + (^e, we can rewrite Eq. 1 1 in terms of Vm and 0e- 

dVm 

Cm V ’ DiS/Vm '^ion ix^m) ~ ^ V (j)e H“ '<^app i^^ H . (12) 

In this form, we have in Ff a system composed of a nonlinear parabolic equation (Eq. 
12) in Vm, coupled with an elliptic equation (Eq. 9) in 0e- 

Under the case of equal anisotropy ratio [10], = kDi , 0e = + 1), 

and Eq. 12 becomes 



dVm 
' dt 



'^ion('Cm) H“ '^app 



(13) 



while in B we have an elliptic equation (Eq. 5) in 0o- 



2.3 Boundary Conditions 

Let the closed surface Sh be a boundary separating bidomain region H and surrounding 
volume conductor B, and let the closed surface Sb bound region B; n will denote the 
unit outward normal to Sh and Sb- Since the potential must be continuous at each 
boundary and the current must be continuous across each boundary, 

0o = 0e, n • (7oV0o = n • (L)V0e + AVi;,n) on Sh (14) 

n -(70X700 = 0 on Sb- (15) 

The bidomain system in H and the passive volume conductor B are connected via 
transmission condition on the boundary Sh (Eq. 14). Since the sources in H are related 
to the presence of intracellular medium, which is absent in B, we may assume that 
the vector DiVvm in Eq. 14 is tangent to the surface Sh', this becomes an additional 
boundary condition: 



n • DiVvm = 0 on Sh- 
The same argument can be made for the vector DiV(j)i 


(16) 


n • DiV(j)i = 0 on Sh- 


(17) 



The conditions 16 and 17 are analogous to the sealed-end condition in the one- 
dimensional cable model. 
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2.4 Extracardiac Fields 

As presented above, the extracardiac potential 0o can be solved everywhere in B ac- 
cording to Eq. 5 and the boundary conditions 14 and 15. The external magnetic field 
B(r) due to the total current density J from Eq. 3 can in turn be evaluated from the 
Biot-Savart law. 



3 Solution of the Bidomain Equations 

3.1 Dalhousie Propagation Model 

Our propagation model was developed at Dalhousie University to simulate the electrical 
excitation wavefronts in two- and three-dimensional arrays of excitable cells [4]. Besides 
the anisotropic bidomain theory, the model comprises features from the cellular automata 
theory. To develop a model for simulating propagated excitation in large and complex 
cardiac network structure, the bisyncytium was tessellated into n subregions (n ^ 10^). 
Each subregion is characterized by two static parameters: a cell type and a principal fiber 
direction a. Associated with each cell is a distinct set of electrophysiological parameters. 
It is assumed that the cell membrane has a voltage threshold Vth- Electrotonic interaction 
of the model’s excitable elements is governed by the parabolic equation 13. 

Formally, the model can be defined as a cellular automaton consisting of n intercon- 
nected finite state machines. The time domain is discretized into a set of time instants, 
t/c+i = tk~\- St, where St is the discrete time step. Each cell is assigned its own space in 
the computer memory, where both static and dynamic information are stored. The static 
information specifies the cell type and the local fiber direction a (in terms of the angles 
0 in the globally defined polar coordinate system). The cell’s dynamic information 
consists of two components: the macro-state and the micro-state. Each cell can be in one 
of the four macro-states, which correspond to well-known distinct phases of the action 
potential: 1) the resting state, 2) the excitatory state, 3) the absolute refractory state, 4) 
the relative refractory state. There are 16 possible macro-state transitions between the 
four states, but physiological constraints reduce possible transitions to eleven. The states 
and the transitions are described in detail elsewhere [4,5]. 



3.2 Equal Anisotropy Ratio 

The membrane potential (t + 1 ) is solved from discretized Eq. 1 3 by explicit difference 

method for Vm{t) [4,13]. Presently, the ionic current density iion in Eq. 13 includes only 
the inward rectifier ixi [14]. The suprathreshold behavior of the elements depends only 
on time, the recovery interval and the cell type; it is largely determined by the state- 
transition rules of the cellular automaton [4,5]. 

The model has been previously tested in simulating the propagation in two- and 
three-dimensional arrays of cells [6,12], and also in a realistic geometry model of the 
human ventricles for simulating potential maps and magnetic field distributions of normal 
ventricular depolarization (Horacek et ah, this volume). 
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3.3 General Anisotropy Ratio 

In a general case, we have to solve the parabolic Eq. 12 and the elliptic Eq. 9, with 
the appropriate boundary conditions 14 and 16. The solution becomes computationally 
tedious, but can be achieved with utilizing the sparse nature of the discretized problem. 
For example, Roth [15] and Henriquez et al. [16] solved the elliptic equation 9 using an 
iterative method, and then substituted the resulting (j)^ to Eq. 12. 

The bidomain equations are solvable also in the context of Dalhousie hybrid model. 
If we rewrite Eq. 4 as V • {DiVv^ + DV(j)e) = 0, Eq. 12 becomes 

Ov 

Qf ~ ^app '^ion('^m) ^ 5 (1^) 

which needs to be solved simultaneously with Eq. 9: 

V • DV(l)e = -V • DiVvm . (19) 

1. Initially, Vm{t = 0) = Vr in all cells, and selecting = 0, (j)e{t = 0) = —Vr. 

2. When an external stimulus is applied at time tk, solve as a function of 

4>e{tk), applying the explicit forward method to parabolic Eq. 18. 

3. Substitute the resulting to the elliptic Eq. 19, and solve the distribution of 

0e(^/c+l)- 

4. Switch to the next time point and repeat the previous two steps. 

The fortran-language implementation was run on an IBM RS/6000 SP supercomputer at 
the CSC - Scientific Computing Ltd (Espoo, Finland; http://www.csc.fi). In solving the 
parabolic Eq. 18, we utilized the parallelized sparse system solvers of the ESSE library 
(http://www.rs6000.ibm.com/resource/aix_resource/sp_books/essl/). 

4 Simulations 

We performed computer simulations in an array of 100 x 100 x 10 cells with the spacing 
ofh = 0.2 mm. Initial stimulus was delivered at the center of the array and the excitation 
was propagated for 20 ms. Figure 1 displays the first test simulations for excitation 
isochrones (i.e., the times when Vm at each cell exceeds the threshold to trigger an 
action potential) for two sets of conductivities: 

- For equal anisotropy ratio, Eq. 13, = 2.5 mS/cm, = 0.6 mS/cm. 

- For unequal anisotropy ratio, Eqs. 18 and 19, = 3 mS/cm, aj = 0.315 mS/cm, 

af = 2 mS/cm, = 1.351 mS/cm (according to [16] and [17]). 

- Other parameters: 6t = 0.01 ms. Cm = 1.0 /iF/cm^, x = 2000. 

In the equal anisotropy ratio case, the axial and transversal propagation velocities 
were, respectively, about 1.0 and 0.5 mm/ms. The propagation ellipsoid for the general 
anisotropy showed a slightly ’slimmer’ shape, the axial propagation velocity was about 
the same but the transversal propagation was slightly slower than in the equal anisotropy 



case. 
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Fig. 1. Simulated isochrones for a) equal anisotropy ratio, and for b) unequal anisotropy ratio. The 
simulation parameter values are given in the text. 



5 Discussion 

In this study, the viewpoint is macroscopic, and each excitable element actually rep- 
resents a large number (^1000) of myocardial cells. We used a fixed shape for the 
action potential above the threshold value, but the shape and duration can easily be al- 
tered to match experimentally measured action potentials or action potentials simulated 
with one- or two-dimensional cellular-scale models based on Luo-Rudy-type differential 
equations [14]. 

Previous studies of general anisotropy ratios (e.g. [15], [16]) were involving detailed 
characteristics of the ionic channels in fairly small blocks of excitable cells {h ^ lOO/zm). 
Such approach may not be easily extendable to a macroscopic view required in simula- 
tions of body surface maps. Our hybrid approach can in turn accommodate different grid 
sizes when special care is paid in defining the model parameters and difference formula 
[13]. 

The isochrones in our first simulations with general anisotropy ratios have elliptical 
shapes, and the estimated average propagation velocities are in fairly good agreement 
with the square-root relation [1^49]. Moreover, the velocities are in 

agreement with the velocities estimated from in vitro experimental data [20]. 

Our model allows simulations of propagated excitation in large and complex anatom- 
ical structures, such as the ventricles, by assuming equal anisotropy ratio and bypassing 
accurate modelling of membrane ionic processes during the cardiac action potential. 
More general anisotropies can also be taken into account by increasing the computational 
demands (intercoupled elliptic and parabolic PDE). Still, more detailed simulations are 
required to study the effect of different anisotropy ratios, especially on body surface 
maps, and to validate the assumptions for the whole-heart simulations. 
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Abstract. Activation time (AT) imaging from electrocardiographic 
mapping data is under development becoming feasible for clinical ap- 
plications. In this study the AT imaging approach is applied to data of 
the human atria and ventricles in four patients. The reconstructed atrial 
AT pattern on the endo- and epicardium for paced rhythm data was 
compared with the CARTO® map. The geometrical error between the 
modelled right atrial (RA) target chamber and the anatomical CARTO® 
map was between 4 and 7mm. With regard to the localization error we 
compared the locations of the first (endocardial) breakthroughs deter- 
mined from the reconstructed AT map and from the CARTO® map 
within RA. This localization error was identified to be within 6 and 
13mm. 



1 Introduction 

Combined magnetic resonance imaging (MRI) of the torso anatomy and elec- 
trocardiographic (ECG) mapping enables noninvasive imaging of the electrical 
function in the human heart [1,2, 3, 4, 5]. Beside the imaging of cardiac movement, 
perfusion and metabolism, the reconstruction of electrical function will have a 
significant clinical impact. 

The primary source in the cardiac muscle is the spatiotemporal transmem- 
brane potential iprn distribution. For depolarization, the assumption of electrical 
isotropy in the myocardium yields reasonable results [6,7]. Because of the un- 
known individual fiber orientation this hypothesis is considered in the uniform 
dipole layer theory based source model formulation. Herewith, the forward and 
inverse formulation can be reduced to a two-dimensional field and scalar poten- 
tial problem. In general, the boundary element method is applied for this kind 
of problem [6]. 

In the inverse problem parameters describing features of ^pm or the epicardial 
(more precisely the pericardial) potential are estimated [1,2, 4, 8, 5]. The epicardial 
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potential as well as the potential on all other conductivity interfaces are related 
to iprn by a Fredholm integral equation of 2^^ order [6,7]. The most established 
inverse formulations are the imaging of the activation time (AT) map on the 
entire surface of the ventricle or the atrium (employing the uniform dipole layer 
or transmembrane potential source model formulation) and of the epicardial po- 
tential pattern (by solving the so-called epicardial potential problem) [1,2, 9, 4, 8]. 

In clinical electrocardiology arrhythmias, in particular within the human 
atrium and their therapeutic treatment, play a major role [10]. In clinical prac- 
tise, the localization of the origin of arrhythmias is currently achieved by tra- 
ditional catheter techniques and by recently introduced catheter mapping tech- 
niques, e.g., by CARTO® (Biosense Webster Inc., A Johnsonfe Johnson com- 
pany) [11]. These methods have limitations, e.g., they do not allow acquiring a 
single beat activation map and are to some extent time consuming. The single 
beat AT imaging approach would provide this information fully noninvasive and 
- at least theoretically - immediately after the ECG recording. From a clini- 
cal point of view this might have several advantages. The AT imaging methods 
permit the reconstruction of single focal, multiple focal, and more distributed 
activation patterns. Furthermore, these methods can distinguish between areas 
with early and late activation. Potential clinical applications of the ECG in- 
verse problem are noninvasive imaging of atrial and ventricular ectopic as well 
as pre-excited activation. 

In this study we investigated the applicability of the AT imaging method in 
four patients who underwent an electrophysiological (EP) study in the catheter 
laboratory. We focused on paced atrial and ventricular activation sequences. 

2 Methods 

Prior to the treatment in the catheter laboratory individual anatomical data 
were obtained by MRI using a Magnetom- Vision-Plus 1.5 Tesla scanner. Atrial 
and ventricular geometry was recorded in CINE-mode during breath-hold (expi- 
ration, 21 X 7 oblique short axis scans, 4 and 6mm spacing). The lungs and the 
torso shape were recorded in Tl-ELASH-mode during breath-hold (expiration, 
40 axial scans, 10mm spacing). 12 markers (vitamin E capsules, 7 anatomical 
landmarks on the anterior and lateral chest wall, 5 electrode positions on the 
patient’s back) were used to couple the electrode positions and the reference 
points to the MRI frame. Erom this data set a boundary element volume con- 
ductor model was built up. Such a model is depicted in Eig. 1 from an anterior left 
lateral view for a 21-year old male patient suffering from the Wolff- Parkinson- 
White(WPW)-syndrome. 

The patients were then moved to the catheter laboratory and EGG mapping 
data were recorded during and after the EP study. ECG mapping data were 
collected in 62-channels by the Mark-8 system (Biosemi V.O.E., Amsterdam, The 
Netherlands). A Wilson terminal defined the reference potential. The sampling 
rate was 2048 Hz. Signals were bandpass filtered with a lower and upper edge 
frequency of 0.3 Hz and 400 Hz, respectively. The AC-resolution of the system is 
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500 nV/bit (16 bit per channel). Radiotransparent carbon electrodes were used 
in order to allow simultaneous X-ray examination. 

The position of the electrodes on the anterior and lateral chest wall was 
digitized by the Fastrak® system (Polhemus Inc., Colchester, Vermont). Addi- 
tionally, the positions of the 7 anterior and lateral landmarks were digitized in 
order to allow coordinate transformation to the MRI frame. The locations of the 
5 upper posterior electrodes were identical with the position of the 5 posterior 
MRI markers. 




Fig. 1. Volume conductor model of a 21-year old male patient suffering from the WPW- 
syndrome. From the patient axial MRI scans for modelling the chest and the right and 
left lung, and ECG-gated short axis scans through the cardiac muscle for modelling 
the atrium (600ms) and the ventricle (0ms) were acquired. For this patient 45 elec- 
trodes (spheres) were considered for reconstructing the atrial and ventricular AT map. 
The ventricular and atrial surface model consisted of 616 and 674 nodal points (= 
unknowns), respectively. The Wilson terminal was considered as reference. The elec- 
trode and reference locations were acquired with the Fastrak® system in the catheter 
laboratory. 



The EGG raw data were pre-processed by baseline correction, but no ad- 
ditional filtering was applied [12]. The transfer matrix was calculated applying 
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the boundary element method with linear triangular elements considering the 
Wilson terminal [6]. The AT map was determined from single beat ECG map- 
ping data applying an iterative linear inverse approach. Details on this inverse 
approach and on the AT imaging approach in general, can be found elsewhere 
[1,2,4]. 

Because of the complex anatomy of the human atrium there are several as- 
pects in the evaluation of a ’’good” atrial surface model. Specific care has to be 
taken on the MRI protocol in obtaining scans with a high gray- value contrast. 
The atrial surface model shown in Fig. 2 and 3 was obtained from short-axis 
scans with 4mm slice thickness in a CINE-mode during breath-hold. For the 
atrial end-diastolic phase the trigger delay was 600ms after the R-peak. The 
epicardium was modelled assuming a uniform wall thickness of 4mm for the left 
and right atrial free wall. Specific attention has to be paid to the identification 
of the pulmonary veins, the vena cava inferior and superior as well as to the 
tricuspid and mitral annulus. The modeling of these structures is important for 
this kind of inverse formulation. Because of the sophisticated curvature and of 
the narrowness of the individual endo- and epicardial boundary element surfaces 
the proper meshing is of utmost importance in order to obtain a high-quality 
transfer matrix. 

3 Results 

In addition to the 62-channel ECG mapping and MRI data a CARTO® map 
for the right atrium (RA) was available from each of the four patients. Two of 
these maps were acquired during sinus rhythm, the other ones during a coronary 
sinus (CS) pacing protocol. AT imaging was performed from single beat sinus 
and from CS and high right atrium (HRA) paced rhythm data. The geometrical 
error between the modelled RA target chamber and the anatomical CARTO® 
map was between 4 and 7mm. With regard to the localization error we compared 
the locations of the first (endocardial) breakthroughs determined from the re- 
constructed AT map and from the CARTO® map. This localization error was 
identified to be within 6 and 13mm. 

In order to give an example of the localization accuracy, two different atrial 
AT map reconstructions achieved from the 21-year old male patient’s data set are 
shown in Fig. 2 and Fig 3. AT imaging was performed for HRA and CS pacing. 
Reconstruction was done from single beat data. For this localization protocol 
the ECG mapping data and the corresponding CARTO® maps were acquired 
at the end of the EP study, after successful ablation of the accessory pathway. 

Firstly, a pacing catheter was placed in the HRA of RA, creating a virtual 
triangle with RAA and VCS. The reconstructed AT map is depicted in Fig. 2. 
It can be seen clearly, that the origin of this paced propagation wave is prop- 
erly found. In this case, however, only an anatomical marker was available for 
comparison. Secondly, the pacing catheter was moved to the CS- Ostium location 
and the stimuli protocol was continued. Again, the propagation pattern was re- 
constructed very well from the stimulated target P wave as can be seen in Fig. 3. 
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Fig. 2. Endocardial activation time map in ms for the HRA pacing protocol seen 
from a caudal oblique view. Isochrones are depicted in steps of 10ms. The following 
abbreviations are used: Right atrium (RA), left atrium (LA), right atrial appendage 
(RAA), left atrial appendage (LAA), vena cava inferior (VCI), vena cava superior 
(VCS), tricuspid annulus (TA), mitral annulus (MA) 




Fig. 3. Endocardial activation time map in ms for the CS-Ostium pacing protocol seen 
from a caudal right lateral oblique view. Isochrones are depicted in steps of 10ms. 
Abbreviations as in Eig. 2. 
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For both reconstructions no target wave signal-subtraction was necessary. The 
stimulated P waves were not superimposed by ventricular activity. The relevance 
of this technique for source imaging is discussed in [12]. 

For the CS pacing site the localization error of the first endocardial break- 
through was determined with 13mm by comparing the RA endocardial AT map 
with the CARTO® activation map of the same target chamber. 

In this 21-year old male patient ECG mapping data were also acquired during 
pacing of the apex of the right ventricle (RV). In Fig. 4 the reconstructed AT map 
is depicted on the endocardium (upper panel) and epicardium (lower panel) from 
an anterior posterior view. Again, reconstruction was done from single beat data. 
As can be seen clearly, the activation spreads from the RV pacing site through 
the Purkinje fibre network to the right and left ventricular (LV) free wall. 

In all four patients ventricular AT maps for sinus and WPW-sinus rhythm 
data were reconstructed. The obtained AT maps were found to be in good agree- 
ment with the well-known sinus rhythm in humans and with the clinical findings 
with regard to the localization of the accessory pathways. 

4 Discussion 

We demonstrated that AT imaging within the human atrium and ventricle from 
paced and sinus rhythm ECG mapping data is feasible under clinical condi- 
tions. Of course, the imaging of spontaneous rhythms, like the onset of flutter or 
spontaneous foci from the pulmonary veins, will be another important milestone 
in the development and clinical establishment of this novel diagnostic imag- 
ing technique. To the extent to which this method may be applicable at the 
present development stage and to our experience attained up to now, there are 
some important aspects for a successful imaging of the target activation pat- 
tern. Firstly, not surprisingly, an accurate model of the target source-containing 
surface coupled as exactly as possible in time and space to the ECG data is of 
utmost importance. Secondly, the target ECG wave for which imaging has to be 
performed should clearly represent the underlying activation process. In clinical 
situations this often is not the case, therefore signal-subtraction techniques will 
have a huge methodical impact. Thirdly, the applied inverse approach must have 
the numerical performance to extract a unique and stable inverse solution in the 
presence of a noisy environment. 

In summary, the presented results achieved on atrial and ventricular AT 
imaging are very promising and raise hope that further research will bring this 
new diagnostic tool closer to clinical applications. 
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RVOT 



Fig. 4. Endo- (upper panel) and epicardial (lower panel) activation time map in ms for 
the RV pacing protocol seen from an anterior posterior view. Isochrones are depicted in 
steps of 10ms. The following abbreviations are used: Right ventricle (RV), left ventricle 
(LV), right ventricular outflow tract (RVOT) 
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